From 9ad4db050e3ee06b2a16936ed7baa9856c1a7981 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 26 Mar 2024 20:49:01 -0600 Subject: [PATCH 01/12] start looking into romberg integration --- src/DuctAPE.jl | 2 + src/precomputation/gausskronrod_integrals.jl | 196 ++++++++++++ src/precomputation/integrands.jl | 315 +++++++++++++++++++ src/precomputation/romberg_integrals.jl | 288 +++++++++++++++++ src/types.jl | 25 ++ 5 files changed, 826 insertions(+) create mode 100644 src/precomputation/gausskronrod_integrals.jl create mode 100644 src/precomputation/integrands.jl create mode 100644 src/precomputation/romberg_integrals.jl create mode 100644 src/types.jl diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index 0bdf8fb6..51c8361f 100644 --- a/src/DuctAPE.jl +++ b/src/DuctAPE.jl @@ -21,6 +21,8 @@ using LinearAlgebra: mul!, ldiv!, lu!, NoPivot, issuccess#, factorize # used in using Printf # used when verbose option is selected +using Primes # used in Romberg integration settings + #---------------------------------# # INCLUDES # #---------------------------------# diff --git a/src/precomputation/gausskronrod_integrals.jl b/src/precomputation/gausskronrod_integrals.jl new file mode 100644 index 00000000..52da4960 --- /dev/null +++ b/src/precomputation/gausskronrod_integrals.jl @@ -0,0 +1,196 @@ +#---------------------------------# +# VORTEX # +#---------------------------------# + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function nominal_vortex_panel_integration( + integration_options::GaussKronrod, + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=[0.0; 1.0], + debug=false, +) + + # Define function to integrate + function fsample(t) + return nominal_vortex_induced_velocity_sample( + t, + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=nondimrange, + ) + end + + V, err = quadgk( + fsample, + 0.0, + 1.0; + order=integration_options.order, + maxevals=integration_options.maxevals, + atol=integration_options.atol, + ) + + if debug + return reshape(V, (2, 2)), err + # return V, err + else + return reshape(V, (2, 2)) + # return V + end +end + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function self_vortex_panel_integration( + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=[0.0; 1.0], + debug=false, +) + + # Define function to integrate + function fsample(t) + return self_vortex_induced_velocity_sample( + t, + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=nondimrange, + ) + end + + V, err = quadgk( + fsample, + 0.0, + 0.5, + 1.0; + order=integration_options.order, + maxevals=integration_options.maxevals, + atol=integration_options.atol, + ) + + cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( + controlpoint[2], influence_length + ) + + V .*= influence_length + V[1:2] .+= cache_vec[1] / 2.0 + + if debug + return reshape(V, (2, 2)), err + else + return reshape(V, (2, 2)) + end +end + +#---------------------------------# +# SOURCE # +#---------------------------------# + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function nominal_source_panel_integration( + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=[0.0; 1.0], + debug=false, +) + + # Define function to integrate + function fsample(t) + return nominal_source_induced_velocity_sample( + t, + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=nondimrange, + ) + end + + V, err = quadgk( + fsample, + 0.0, + 1.0; + order=integration_options.order, + maxevals=integration_options.maxevals, + atol=integration_options.atol, + ) + + if debug + return reshape(V, (2, 2)), err + # return V, err + else + return reshape(V, (2, 2)) + # return V + end +end + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function self_source_panel_integration( + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=[0.0; 1.0], + debug=false, +) + + # Define function to integrate + function fsample(t) + return self_source_induced_velocity_sample( + t, + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=nondimrange, + ) + end + + V, err = quadgk( + fsample, + 0.0, + 0.5, + 1.0; + order=integration_options.order, + maxevals=integration_options.maxevals, + atol=integration_options.atol, + ) + + cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( + controlpoint[2], influence_length + ) + + V .*= influence_length + V[3:4] .+= cache_vec[2] / 2.0 + + if debug + return reshape(V, (2, 2)), err + else + return reshape(V, (2, 2)) + end +end diff --git a/src/precomputation/integrands.jl b/src/precomputation/integrands.jl new file mode 100644 index 00000000..820e2b9d --- /dev/null +++ b/src/precomputation/integrands.jl @@ -0,0 +1,315 @@ +###################################################################### +# # +# VORTEX # +# # +###################################################################### + +#---------------------------------# +# Nominal # +#---------------------------------# + +""" + +# Arguments: +- `t::Float` : sample point in range (0,1) selected by quadrature. +""" +function nominal_vortex_induced_velocity_sample( + t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] +) + + # Transform from (0,1) to actual position on panel + cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) # z coordinate + cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) # r coordinate + + # get relative geometry: xi, rho, m, r0 = calculate_xrm(controlpoint, [z; r]) + calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) + + # Get velocity components at sample points + cache_vec[7] = vortex_ring_vz!( + cache_vec[3],#xi + cache_vec[4],#rho + cache_vec[5],#m + cache_vec[6],#rj + 1.0, + view(cache_vec, 11:15), + ) #vz + + cache_vec[8] = vortex_ring_vr!( + cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[6], view(cache_vec, 11:16) + ) #vr + + #= + assemble output components in the format: + [x_j r_j; x_{j+1} r_{j+1}] + and scale by influence panel length + (due to transformation of integration range to/from range=(0,1)) + =# + # return StaticArrays.SMatrix{2,2}( + # [vz * [1.0 - t; t] vr * [1.0 - t; t]] * influence_length + # ) + return StaticArrays.SVector{4}( + [ + cache_vec[7] * (1.0 - t) + cache_vec[7] * t + cache_vec[8] * (1.0 - t) + cache_vec[8] * t + ] * influence_length, + ) + # return [vz * [1.0 - t; t] vr * [1.0 - t; t]] * influence_length +end + +#---------------------------------# +# Singular # +#---------------------------------# + +""" +""" +function subtracted_singular_vortex_influence(node, controlpoint) + rmag2 = (controlpoint[1] - node[1])^2 + (controlpoint[2] - node[2])^2 + num1z = controlpoint[1] - node[1] + num1r = node[2] - controlpoint[2] + den1 = 2.0 * pi * rmag2 + den2 = 64.0 * controlpoint[2]^2 + + if isapprox(rmag2, 0.0) || isapprox(controlpoint[2], 0.0) + axial = 0.0 + else + axial = num1r / den1 - log(rmag2 / den2) / (8.0 * pi * controlpoint[2]) + end + + if isapprox(controlpoint[2], 0.0) + radial = 0.0 + else + radial = num1z / den1 + end + + return axial, radial +end + +function subtracted_singular_vortex_influence!(node, controlpoint, cache_vec) + cache_vec[1] = (controlpoint[1] - node[1])^2 + (controlpoint[2] - node[2])^2 + cache_vec[2] = controlpoint[1] - node[1] + cache_vec[3] = node[2] - controlpoint[2] + cache_vec[4] = 2.0 * pi * cache_vec[1] + cache_vec[5] = 64.0 * controlpoint[2]^2 + + if isapprox(cache_vec[1], 0.0) || isapprox(controlpoint[2], 0.0) + axial = 0.0 + else + axial = + cache_vec[3] / cache_vec[4] - + log(cache_vec[1] / cache_vec[5]) / (8.0 * pi * controlpoint[2]) + end + + if isapprox(controlpoint[2], 0.0) + radial = 0.0 + else + radial = cache_vec[2] / cache_vec[4] + end + + return axial, radial +end + +""" +""" +function analytically_integrated_vortex_influence(r, influence_length) + axial = (influence_length / (4.0 * pi * r)) * (1.0 + log(16.0 * r / influence_length)) + radial = 0.0 + return axial, radial +end + +""" + +# Arguments: +- `t::Float` : sample point in range (0,1) selected by quadrature. +""" +function self_vortex_induced_velocity_sample( + t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] +) + + # Transform from (0,1) to actual position on panel + cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) + cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) + + # get relative geometry + calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) + + # Get velocity components at sample points + cache_vec[7] = vortex_ring_vz!( + cache_vec[3],#xi + cache_vec[4],#rho + cache_vec[5],#m + cache_vec[2],#rj + 1.0, + view(cache_vec, 11:15), + ) #vz + cache_vec[8] = vortex_ring_vr!( + cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[2], view(cache_vec, 11:16) + ) #vr + + # Get singular piece to subtract + cache_vec[9], cache_vec[10] = subtracted_singular_vortex_influence!( + (cache_vec[1], cache_vec[2]), controlpoint, view(cache_vec, 11:15) + ) + + #= + assemble output components in the format: + [x_j r_j; x_{j+1} r_{j+1}] + and scale by influence panel length + (due to transformation of integration range to/from range=(0,1)) + =# + return [ + cache_vec[7] * (1.0 - t) - cache_vec[9] / 2.0 + cache_vec[7] * t .- cache_vec[9] / 2.0 + cache_vec[8] * (1.0 - t) .- cache_vec[10] / 2.0 + cache_vec[8] * t .- cache_vec[10] / 2.0 + ] +end + +###################################################################### +# # +# SOURCE # +# # +###################################################################### + +#---------------------------------# +# Nominal # +#---------------------------------# + +""" + +# Arguments: +- `t::Float` : sample point in range (0,1) selected by quadrature. +""" +function nominal_source_induced_velocity_sample( + t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] +) + + # Transform from (0,1) to actual position on panel + cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) + cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) + + # get relative geometry + calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) + + # Get velocity components at sample points + cache_vec[7] = source_ring_vz!( + cache_vec[3],#xi + cache_vec[4],#rho + cache_vec[5],#m + cache_vec[6],#rj + view(cache_vec, 11:15), + ) + cache_vec[8] = source_ring_vr!( + cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[6], view(cache_vec, 11:16) + ) + + #= + assemble output components in the format: + [x_j r_j; x_{j+1} r_{j+1}] + and scale by influence panel length + (due to transformation of integration range to/from range=(0,1)) + =# + # return StaticArrays.SMatrix{2,2}([cache_vec[7] * [1.0 - t; t] cache_vec[8] * [1.0 - t; t]] * influence_length) + return StaticArrays.SVector{4}( + [ + cache_vec[7] * (1.0 - t) + cache_vec[7] * t + cache_vec[8] * (1.0 - t) + cache_vec[8] * t + ] * influence_length, + ) +end + +#---------------------------------# +# Singular # +#---------------------------------# + +""" +""" +function subtracted_singular_source_influence(node, controlpoint) + rmag2 = (controlpoint[1] - node[1])^2 + (controlpoint[2] - node[2])^2 + #TODO: write up math and check signs here. + num1z = controlpoint[1] - node[1] + num1r = controlpoint[2] - node[2] + den1 = 2.0 * pi * rmag2 + # den2 = 64.0 * controlpoint[2]^2 + den2 = controlpoint[2]^2 #DFDC has no 64 here somehow + + radial = num1r / den1 - log(rmag2 / den2) / (8.0 * pi * controlpoint[2]) + axial = num1z / den1 + + return axial, radial +end + +function subtracted_singular_source_influence!(node, controlpoint, cache_vec) + cache_vec[1] = (controlpoint[1] - node[1])^2 + (controlpoint[2] - node[2])^2 + #TODO: write up math and check signs here. + cache_vec[2] = controlpoint[1] - node[1] + cache_vec[3] = controlpoint[2] - node[2] + cache_vec[4] = 2.0 * pi * cache_vec[1] + # den2 = 64.0 * controlpoint[2]^2 + cache_vec[5] = controlpoint[2]^2 #DFDC has no 64 here somehow + + radial = + cache_vec[3] / cache_vec[4] - + log(cache_vec[1] / cache_vec[5]) / (8.0 * pi * controlpoint[2]) + axial = cache_vec[2] / cache_vec[4] + + return axial, radial +end + +""" +""" +function analytically_integrated_source_influence(r, influence_length) + radial = (influence_length / (4.0 * pi * r)) * (1.0 + log(2.0 * r / influence_length)) + axial = 0.0 + return axial, radial +end + +""" + +# Arguments: +- `t::Float` : sample point in range (0,1) selected by quadrature. +""" +function self_source_induced_velocity_sample( + t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] +) + + # Transform from (0,1) to actual position on panel + cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) + cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) + + # get relative geometry + calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) + + # Get velocity components at sample points + cache_vec[7] = source_ring_vz!( + cache_vec[3],#xi + cache_vec[4],#rho + cache_vec[5],#m + cache_vec[2],#rj + view(cache_vec, 11:15), + ) + cache_vec[8] = source_ring_vr!( + cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[2], view(cache_vec, 11:16) + ) + + # Get singular piece to subtract + cache_vec[9], cache_vec[10] = subtracted_singular_source_influence!( + (cache_vec[1], cache_vec[2]), controlpoint, view(cache_vec, 11:15) + ) + + #= + assemble output components in the format: + [x_j r_j; x_{j+1} r_{j+1}] + and scale by influence panel length + (due to transformation of integration range to/from range=(0,1)) + =# + return [ + cache_vec[7] * (1.0 - t) .- cache_vec[9] / 2.0 + cache_vec[7] * t .- cache_vec[9] / 2.0 + cache_vec[8] * (1.0 - t) .- cache_vec[10] / 2.0 + cache_vec[8] * t .- cache_vec[10] / 2.0 + ] +end diff --git a/src/precomputation/romberg_integrals.jl b/src/precomputation/romberg_integrals.jl new file mode 100644 index 00000000..05a5f29c --- /dev/null +++ b/src/precomputation/romberg_integrals.jl @@ -0,0 +1,288 @@ +#TODO: decide if dfdc method or julia method is better +# the DFDC method requires the whole extrapolation, but you can limit the number of intervals to evaluate when the integral is easy. +# when the integral is hard though, you end up doing 2x more calculations than you would have if you set your sample size larger to begin with. +# the julia method requires you to evaluate all the sample points up front, but will terminate the extrapolation as soon as you hit your desired tolerance. So you need to have a large number of samples to start with. +# DFDC doesn't actually do a true trapezoidal integration, rather it assumes a linear function and does δx*f(midpoint of iterval), this is technically less accurate than the trapezoidal method, but does avoid sampling at the end and midpoints of panels. Theoretically, the extrapolation will account for the less accurate integration method. +# perhaps a combo, where you do a DFDC-like integration (not quite trapezoidal), but do the large sample size and early extrapolation termination. + + +""" + extrapolate!(fh::AbstractVector; power=1, atol=0, rtol=0, maxeval=typemax(Int), breaktol=Inf) + +Similar to `extrapolate(fh)`, performs Richardson extrapolation on an array `fh` +of `(f(h), h)` tuples (in order of decreasing `|h|`), but overwrites the array +`fh` in-place with intermediate calculations. + +(Thus, the array `fh` must be a vector of `Tuple{T,H}` values, where `H<:Number` is +the type of `h` and `T` is the type of the extrapolated `f(0)` **result**. This `T` +should be a floating-point type, i.e. `fh` should contain `float(f(h))` if the +function you are extrapolating is not already floating-point-valued.) +""" +function extrapolate!( + fh::AbstractVector{<:Tuple{TF1,TF2}}; + power::Int=2, + atol::Float64=1e-6, + rtol::Float64=0.0, + breaktol::Float64=Inf, + maxeval::Int=typemax(Int), +) where {TF1,TF2} + + # - Check for Problems - # + (rtol >= 0 && atol >= zero(atol)) || + throw(ArgumentError("rtol and atol must be nonnegative")) + breaktol > 0 || throw(ArgumentError("breaktol must be positive")) + isempty(fh) && throw(ArgumentError("(f,h) array must be non-empty")) + + # - rename for convenience - # + (f0, h) = first(fh) + + # - initialize - # + err::typeof(float(norm(f0))) = Inf + numeval = 1 + maxeval = min(maxeval, length(fh)) + + # - Loop - # + while numeval < maxeval + # - increment - # + numeval += 1 + + # - update - # + (f_prime, h_prime) = fh[numeval + (firstindex(fh) - 1)] + abs(h) > abs(h_prime) || + throw(ArgumentError("|$h_prime| >= |$h| is not decreasing")) + h = h_prime + minerr_prime = oftype(err, Inf) + f_ip1 = f_prime + + # - Inner Loop - # + for i in (numeval - 1):-1:1 + f_i, h_i = fh[i + (firstindex(fh) - 1)] + c = (h_i / h_prime)^power + f_ip1 += (f_ip1 - f_i) / (c - 1) + fh[i + (firstindex(fh) - 1)] = (f_ip1, h_i) + err_prime = norm(f_ip1 - f_i) + minerr_prime = min(minerr_prime, err_prime) + if err_prime < err + f0, err = f_ip1, err_prime + end + end + + # - stop early if error increases too much - # + (minerr_prime > breaktol * err || !isfinite(minerr_prime)) && break + + # - converged - # + err <= max(rtol * norm(f0), atol) && break + end + return (f0, err) +end + +@views function romberg!(fh, step_size, samples, endsum, factors, numfactors; kwargs...) + + # - Do Trapezoidal Integration - # + b, e = firstindex(samples), lastindex(samples) + i = numfactors + 1 + fh[i] = (step_size * (sum(samples[(b + 1):(e - 1)]) + endsum), step_size) + nstep = 1 + for (j, K) in factors + @inbounds for k in 1:K + nstep *= j + sdx = nstep * step_size + if i == 2 # last iteration (empty sum) + fh[1] = (sdx * endsum, sdx) + else + fh[i -= 1] = ( + sdx * (sum(samples[(b + nstep):nstep:(e - nstep)]) + endsum), sdx + ) + end + end + end + + # - Do Richardson Extrapolation until converged - # + return extrapolate!(fh; power=2, kwargs...) +end + +#---------------------------------# +# VORTEX # +#---------------------------------# + +function allocate_integration_containers(dispatch_type::TF, nsamples, nfactors) where {TF} + return (; + err=zeros(TF, 2, 2), #always this size + sample_cache=zeros(TF, 20), #always this size, this is for the intermediate calcs + samples=zeros(TF, nsamples, 4), # always 4 wide + fh=[(TF(1.0), TF(1.0)) for i in 1:(nfactors + 1)], + ) +end + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function nominal_vortex_panel_integration!( + integration_options::Romberg, + V, + node1, + node2, + influence_length, + controlpoint, + containers; + debug=false, +) + + # Define function to integrate + for (s, t) in zip(eachrow(containers.samples), integration_options.samplepoints) + s .= nominal_vortex_induced_velocity_sample( + t, + node1, + node2, + influence_length, + controlpoint, + containers.sample_cache; + nondimrange=integration_options.nondimrange, + ) + end + + for (i, s) in enumerate(eachcol(containers.samples)) + V[i], containers.err[i] = romberg!( + containers.fh, + integration_options.stepsize, + s, + (s[1] + s[end]) / 2.0, + factors, + length(factors); + atol=integration_options.atol, + ) + end + + if debug + # return reshape(V, (2, 2)), err + return V, err + else + # return reshape(V, (2, 2)) + return V + end +end + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function self_vortex_panel_integration( + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=[0.0; 1.0], + debug=false, +) + + # Define function to integrate + function fsample(t) + return self_vortex_induced_velocity_sample( + t, + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=nondimrange, + ) + end + + V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) + + cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( + controlpoint[2], influence_length + ) + + V .*= influence_length + V[1:2] .+= cache_vec[1] / 2.0 + + if debug + return reshape(V, (2, 2)), err + else + return reshape(V, (2, 2)) + end +end + +#---------------------------------# +# SOURCE # +#---------------------------------# + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function nominal_source_panel_integration( + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=[0.0; 1.0], + debug=false, +) + + # Define function to integrate + function fsample(t) + return nominal_source_induced_velocity_sample( + t, + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=nondimrange, + ) + end + + V, err = quadgk(fsample, 0.0, 1.0; order=3, maxevals=1e2, atol=1e-6) + + if debug + return reshape(V, (2, 2)), err + # return V, err + else + return reshape(V, (2, 2)) + # return V + end +end + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function self_source_panel_integration( + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=[0.0; 1.0], + debug=false, +) + + # Define function to integrate + function fsample(t) + return self_source_induced_velocity_sample( + t, + node1, + node2, + influence_length, + controlpoint, + cache_vec; + nondimrange=nondimrange, + ) + end + + V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) + + cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( + controlpoint[2], influence_length + ) + + V .*= influence_length + V[3:4] .+= cache_vec[2] / 2.0 + + if debug + return reshape(V, (2, 2)), err + else + return reshape(V, (2, 2)) + end +end diff --git a/src/types.jl b/src/types.jl new file mode 100644 index 00000000..72d13478 --- /dev/null +++ b/src/types.jl @@ -0,0 +1,25 @@ +abstract type IntegrationMethod end + +@kwdef struct Romberg <: IntegrationMethod + nsamples::Int = 1025 + stepsize::Float64 = 1.0 / 1024 + factors::Primes.Factorization{Int} = Primes.factor(1024) + nfactors::Int = 10 + atol::Float64 = 1e-6 +end + +function set_romberg_options(; nsamples=1025, atol=1e-6) + stepsize = 1.0 / (nsamples - 1) + factors = Primes.factor(nsamples - 1) + nfactors = sum(values(factors)) + + return Romberg(; + atol=atol, nsamples=nsamples, factors=factors, nfactors=nfactors, stepsize=stepsize + ) +end + +@kwdef struct GaussKronrod <: IntegrationMethod + order::Int = 3 + maxevals::Int = 1e2 + atol::Float64 = 1e-6 +end From c88fe972b7d25a412a56659e7a07781a3fa11b7b Mon Sep 17 00:00:00 2001 From: juddmehr Date: Wed, 27 Mar 2024 14:30:38 -0600 Subject: [PATCH 02/12] fix types to be concrete and reduce allocations --- Project.toml | 1 + src/DuctAPE.jl | 3 + src/precomputation/integrals.jl | 296 ------------------------ src/precomputation/integrands.jl | 42 ++++ src/precomputation/romberg_integrals.jl | 191 ++++++++------- src/types.jl | 21 +- 6 files changed, 150 insertions(+), 404 deletions(-) diff --git a/Project.toml b/Project.toml index 4c964776..674fa3ab 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index 51c8361f..4a9bebbc 100644 --- a/src/DuctAPE.jl +++ b/src/DuctAPE.jl @@ -30,6 +30,7 @@ using Primes # used in Romberg integration settings ##### ----- UTILITIES ----- ##### # general utility functions include("utilities/utils.jl") +include("types.jl") ##### ----- PRECOMPUTATION ----- ##### # Pre-solve initializations @@ -52,6 +53,8 @@ include("precomputation/velocities.jl") include("precomputation/body_aic.jl") include("precomputation/induced_velocity_matrices.jl") +include("precomputation/romberg_integrals.jl") +include("precomputation/integrands.jl") ##### ----- SOLVER ----- ##### include("solve/solve.jl") diff --git a/src/precomputation/integrals.jl b/src/precomputation/integrals.jl index ec34c50f..080b0204 100644 --- a/src/precomputation/integrals.jl +++ b/src/precomputation/integrals.jl @@ -8,56 +8,6 @@ # Nominal # #---------------------------------# -""" - -# Arguments: -- `t::Float` : sample point in range (0,1) selected by quadrature. -""" -function nominal_vortex_induced_velocity_sample( - t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] -) - - # Transform from (0,1) to actual position on panel - cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) # z coordinate - cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) # r coordinate - - # get relative geometry: xi, rho, m, r0 = calculate_xrm(controlpoint, [z; r]) - calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) - - # Get velocity components at sample points - cache_vec[7] = vortex_ring_vz!( - cache_vec[3],#xi - cache_vec[4],#rho - cache_vec[5],#m - cache_vec[6],#rj - 1.0, - view(cache_vec, 11:15), - ) #vz - - cache_vec[8] = vortex_ring_vr!( - cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[6], view(cache_vec, 11:16) - ) #vr - - #= - assemble output components in the format: - [x_j r_j; x_{j+1} r_{j+1}] - and scale by influence panel length - (due to transformation of integration range to/from range=(0,1)) - =# - # return StaticArrays.SMatrix{2,2}( - # [vz * [1.0 - t; t] vr * [1.0 - t; t]] * influence_length - # ) - return StaticArrays.SVector{4}( - [ - cache_vec[7] * (1.0 - t) - cache_vec[7] * t - cache_vec[8] * (1.0 - t) - cache_vec[8] * t - ] * influence_length, - ) - # return [vz * [1.0 - t; t] vr * [1.0 - t; t]] * influence_length -end - """ `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ @@ -84,10 +34,8 @@ function nominal_vortex_panel_integration( ) end - V, err = quadgk(fsample, 0.0, 1.0; order=3, maxevals=1e2, atol=1e-6) - if debug return reshape(V, (2, 2)), err # return V, err @@ -101,110 +49,6 @@ end # Singular # #---------------------------------# -""" -""" -function subtracted_singular_vortex_influence(node, controlpoint) - rmag2 = (controlpoint[1] - node[1])^2 + (controlpoint[2] - node[2])^2 - num1z = controlpoint[1] - node[1] - num1r = node[2] - controlpoint[2] - den1 = 2.0 * pi * rmag2 - den2 = 64.0 * controlpoint[2]^2 - - if isapprox(rmag2, 0.0) || isapprox(controlpoint[2], 0.0) - axial = 0.0 - else - axial = num1r / den1 - log(rmag2 / den2) / (8.0 * pi * controlpoint[2]) - end - - if isapprox(controlpoint[2], 0.0) - radial = 0.0 - else - radial = num1z / den1 - end - - return axial, radial -end - -function subtracted_singular_vortex_influence!(node, controlpoint, cache_vec) - cache_vec[1] = (controlpoint[1] - node[1])^2 + (controlpoint[2] - node[2])^2 - cache_vec[2] = controlpoint[1] - node[1] - cache_vec[3] = node[2] - controlpoint[2] - cache_vec[4] = 2.0 * pi * cache_vec[1] - cache_vec[5] = 64.0 * controlpoint[2]^2 - - if isapprox(cache_vec[1], 0.0) || isapprox(controlpoint[2], 0.0) - axial = 0.0 - else - axial = - cache_vec[3] / cache_vec[4] - - log(cache_vec[1] / cache_vec[5]) / (8.0 * pi * controlpoint[2]) - end - - if isapprox(controlpoint[2], 0.0) - radial = 0.0 - else - radial = cache_vec[2] / cache_vec[4] - end - - return axial, radial -end - -""" -""" -function analytically_integrated_vortex_influence(r, influence_length) - axial = (influence_length / (4.0 * pi * r)) * (1.0 + log(16.0 * r / influence_length)) - radial = 0.0 - return axial, radial -end - -""" - -# Arguments: -- `t::Float` : sample point in range (0,1) selected by quadrature. -""" -function self_vortex_induced_velocity_sample( - t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] -) - - # Transform from (0,1) to actual position on panel - cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) - cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) - - # get relative geometry - calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) - - # Get velocity components at sample points - cache_vec[7] = vortex_ring_vz!( - cache_vec[3],#xi - cache_vec[4],#rho - cache_vec[5],#m - cache_vec[2],#rj - 1.0, - view(cache_vec, 11:15), - ) #vz - cache_vec[8] = vortex_ring_vr!( - cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[2], view(cache_vec, 11:16) - ) #vr - - # Get singular piece to subtract - cache_vec[9], cache_vec[10] = subtracted_singular_vortex_influence!( - (cache_vec[1], cache_vec[2]), controlpoint, view(cache_vec, 11:15) - ) - - #= - assemble output components in the format: - [x_j r_j; x_{j+1} r_{j+1}] - and scale by influence panel length - (due to transformation of integration range to/from range=(0,1)) - =# - return [ - cache_vec[7] * (1.0 - t) - cache_vec[9] / 2.0 - cache_vec[7] * t .- cache_vec[9] / 2.0 - cache_vec[8] * (1.0 - t) .- cache_vec[10] / 2.0 - cache_vec[8] * t .- cache_vec[10] / 2.0 - ] -end - """ `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ @@ -231,10 +75,8 @@ function self_vortex_panel_integration( ) end - V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) - cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( controlpoint[2], influence_length ) @@ -259,51 +101,6 @@ end # Nominal # #---------------------------------# -""" - -# Arguments: -- `t::Float` : sample point in range (0,1) selected by quadrature. -""" -function nominal_source_induced_velocity_sample( - t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] -) - - # Transform from (0,1) to actual position on panel - cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) - cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) - - # get relative geometry - calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) - - # Get velocity components at sample points - cache_vec[7] = source_ring_vz!( - cache_vec[3],#xi - cache_vec[4],#rho - cache_vec[5],#m - cache_vec[6],#rj - view(cache_vec, 11:15), - ) - cache_vec[8] = source_ring_vr!( - cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[6], view(cache_vec, 11:16) - ) - - #= - assemble output components in the format: - [x_j r_j; x_{j+1} r_{j+1}] - and scale by influence panel length - (due to transformation of integration range to/from range=(0,1)) - =# - # return StaticArrays.SMatrix{2,2}([cache_vec[7] * [1.0 - t; t] cache_vec[8] * [1.0 - t; t]] * influence_length) - return StaticArrays.SVector{4}( - [ - cache_vec[7] * (1.0 - t) - cache_vec[7] * t - cache_vec[8] * (1.0 - t) - cache_vec[8] * t - ] * influence_length, - ) -end - """ `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ @@ -330,10 +127,8 @@ function nominal_source_panel_integration( ) end - V, err = quadgk(fsample, 0.0, 1.0; order=3, maxevals=1e2, atol=1e-6) - if debug return reshape(V, (2, 2)), err # return V, err @@ -347,95 +142,6 @@ end # Singular # #---------------------------------# -""" -""" -function subtracted_singular_source_influence(node, controlpoint) - rmag2 = (controlpoint[1] - node[1])^2 + (controlpoint[2] - node[2])^2 - #TODO: write up math and check signs here. - num1z = controlpoint[1] - node[1] - num1r = controlpoint[2] - node[2] - den1 = 2.0 * pi * rmag2 - # den2 = 64.0 * controlpoint[2]^2 - den2 = controlpoint[2]^2 #DFDC has no 64 here somehow - - radial = num1r / den1 - log(rmag2 / den2) / (8.0 * pi * controlpoint[2]) - axial = num1z / den1 - - return axial, radial -end - -function subtracted_singular_source_influence!(node, controlpoint, cache_vec) - cache_vec[1] = (controlpoint[1] - node[1])^2 + (controlpoint[2] - node[2])^2 - #TODO: write up math and check signs here. - cache_vec[2] = controlpoint[1] - node[1] - cache_vec[3] = controlpoint[2] - node[2] - cache_vec[4] = 2.0 * pi * cache_vec[1] - # den2 = 64.0 * controlpoint[2]^2 - cache_vec[5] = controlpoint[2]^2 #DFDC has no 64 here somehow - - radial = - cache_vec[3] / cache_vec[4] - - log(cache_vec[1] / cache_vec[5]) / (8.0 * pi * controlpoint[2]) - axial = cache_vec[2] / cache_vec[4] - - return axial, radial -end - -""" -""" -function analytically_integrated_source_influence(r, influence_length) - radial = (influence_length / (4.0 * pi * r)) * (1.0 + log(2.0 * r / influence_length)) - axial = 0.0 - return axial, radial -end - -""" - -# Arguments: -- `t::Float` : sample point in range (0,1) selected by quadrature. -""" -function self_source_induced_velocity_sample( - t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] -) - - # Transform from (0,1) to actual position on panel - cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) - cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) - - # get relative geometry - calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) - - # Get velocity components at sample points - cache_vec[7] = source_ring_vz!( - cache_vec[3],#xi - cache_vec[4],#rho - cache_vec[5],#m - cache_vec[2],#rj - view(cache_vec, 11:15), - ) - cache_vec[8] = source_ring_vr!( - cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[2], view(cache_vec, 11:16) - ) - - # Get singular piece to subtract - cache_vec[9], cache_vec[10] = subtracted_singular_source_influence!( - (cache_vec[1], cache_vec[2]), controlpoint, view(cache_vec, 11:15) - ) - - #= - assemble output components in the format: - [x_j r_j; x_{j+1} r_{j+1}] - and scale by influence panel length - (due to transformation of integration range to/from range=(0,1)) - =# - return [ - cache_vec[7] * (1.0 - t) .- cache_vec[9] / 2.0 - cache_vec[7] * t .- cache_vec[9] / 2.0 - cache_vec[8] * (1.0 - t) .- cache_vec[10] / 2.0 - cache_vec[8] * t .- cache_vec[10] / 2.0 - ] -end - """ `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ @@ -462,10 +168,8 @@ function self_source_panel_integration( ) end - V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) - cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( controlpoint[2], influence_length ) diff --git a/src/precomputation/integrands.jl b/src/precomputation/integrands.jl index 820e2b9d..65275aff 100644 --- a/src/precomputation/integrands.jl +++ b/src/precomputation/integrands.jl @@ -7,6 +7,48 @@ #---------------------------------# # Nominal # #---------------------------------# +""" + +# Arguments: +- `t::Float` : sample point in range (0,1) selected by quadrature. +""" +function nominal_vortex_induced_velocity_sample!( + V, t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] +) + + # Transform from (0,1) to actual position on panel + cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) # z coordinate + cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) # r coordinate + + # get relative geometry: xi, rho, m, r0 = calculate_xrm(controlpoint, [z; r]) + calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) + + # Get velocity components at sample points + cache_vec[7] = vortex_ring_vz!( + cache_vec[3],#xi + cache_vec[4],#rho + cache_vec[5],#m + cache_vec[6],#rj + 1.0, + view(cache_vec, 11:15), + ) #vz + + cache_vec[8] = vortex_ring_vr!( + cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[6], view(cache_vec, 11:16) + ) #vr + + #= + assemble output components in the format: + [x_j r_j; x_{j+1} r_{j+1}] + and scale by influence panel length + (due to transformation of integration range to/from range=(0,1)) + =# + V[1] = cache_vec[7] * (1.0 - t) * influence_length + V[2] = cache_vec[7] * t * influence_length + V[3] = cache_vec[8] * (1.0 - t) * influence_length + V[4] = cache_vec[8] * t * influence_length + return V +end """ diff --git a/src/precomputation/romberg_integrals.jl b/src/precomputation/romberg_integrals.jl index 05a5f29c..214ea9d0 100644 --- a/src/precomputation/romberg_integrals.jl +++ b/src/precomputation/romberg_integrals.jl @@ -5,6 +5,7 @@ # DFDC doesn't actually do a true trapezoidal integration, rather it assumes a linear function and does δx*f(midpoint of iterval), this is technically less accurate than the trapezoidal method, but does avoid sampling at the end and midpoints of panels. Theoretically, the extrapolation will account for the less accurate integration method. # perhaps a combo, where you do a DFDC-like integration (not quite trapezoidal), but do the large sample size and early extrapolation termination. +# Current Plan: do romberg for the non-singular ones, and quadgk for the singular ones. """ extrapolate!(fh::AbstractVector; power=1, atol=0, rtol=0, maxeval=typemax(Int), breaktol=Inf) @@ -78,7 +79,7 @@ end @views function romberg!(fh, step_size, samples, endsum, factors, numfactors; kwargs...) - # - Do Trapezoidal Integration - # + # - Do Trapezoidal Integration with Composite Trapezoidal Rule - # b, e = firstindex(samples), lastindex(samples) i = numfactors + 1 fh[i] = (step_size * (sum(samples[(b + 1):(e - 1)]) + endsum), step_size) @@ -101,10 +102,6 @@ end return extrapolate!(fh; power=2, kwargs...) end -#---------------------------------# -# VORTEX # -#---------------------------------# - function allocate_integration_containers(dispatch_type::TF, nsamples, nfactors) where {TF} return (; err=zeros(TF, 2, 2), #always this size @@ -114,6 +111,10 @@ function allocate_integration_containers(dispatch_type::TF, nsamples, nfactors) ) end +#---------------------------------# +# VORTEX # +#---------------------------------# + """ `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ @@ -130,14 +131,8 @@ function nominal_vortex_panel_integration!( # Define function to integrate for (s, t) in zip(eachrow(containers.samples), integration_options.samplepoints) - s .= nominal_vortex_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - containers.sample_cache; - nondimrange=integration_options.nondimrange, + nominal_vortex_induced_velocity_sample!( + s, t, node1, node2, influence_length, controlpoint, containers.sample_cache ) end @@ -147,8 +142,8 @@ function nominal_vortex_panel_integration!( integration_options.stepsize, s, (s[1] + s[end]) / 2.0, - factors, - length(factors); + integration_options.factors, + integration_options.nfactors; atol=integration_options.atol, ) end @@ -162,47 +157,47 @@ function nominal_vortex_panel_integration!( end end -""" -`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -""" -function self_vortex_panel_integration( - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=[0.0; 1.0], - debug=false, -) - - # Define function to integrate - function fsample(t) - return self_vortex_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, - ) - end - - V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) - - cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( - controlpoint[2], influence_length - ) - - V .*= influence_length - V[1:2] .+= cache_vec[1] / 2.0 - - if debug - return reshape(V, (2, 2)), err - else - return reshape(V, (2, 2)) - end -end +# """ +# `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +# """ +# function self_vortex_panel_integration( +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=[0.0; 1.0], +# debug=false, +# ) + +# # Define function to integrate +# function fsample(t) +# return self_vortex_induced_velocity_sample( +# t, +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=nondimrange, +# ) +# end + +# V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) + +# cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( +# controlpoint[2], influence_length +# ) + +# V .*= influence_length +# V[1:2] .+= cache_vec[1] / 2.0 + +# if debug +# return reshape(V, (2, 2)), err +# else +# return reshape(V, (2, 2)) +# end +# end #---------------------------------# # SOURCE # @@ -211,7 +206,7 @@ end """ `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ -function nominal_source_panel_integration( +function nominal_source_panel_integration!( node1, node2, influence_length, @@ -245,44 +240,44 @@ function nominal_source_panel_integration( end end -""" -`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -""" -function self_source_panel_integration( - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=[0.0; 1.0], - debug=false, -) - - # Define function to integrate - function fsample(t) - return self_source_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, - ) - end - - V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) - - cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( - controlpoint[2], influence_length - ) - - V .*= influence_length - V[3:4] .+= cache_vec[2] / 2.0 - - if debug - return reshape(V, (2, 2)), err - else - return reshape(V, (2, 2)) - end -end +# """ +# `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +# """ +# function self_source_panel_integration( +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=[0.0; 1.0], +# debug=false, +# ) + +# # Define function to integrate +# function fsample(t) +# return self_source_induced_velocity_sample( +# t, +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=nondimrange, +# ) +# end + +# V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) + +# cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( +# controlpoint[2], influence_length +# ) + +# V .*= influence_length +# V[3:4] .+= cache_vec[2] / 2.0 + +# if debug +# return reshape(V, (2, 2)), err +# else +# return reshape(V, (2, 2)) +# end +# end diff --git a/src/types.jl b/src/types.jl index 72d13478..ba59725f 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,11 +1,12 @@ abstract type IntegrationMethod end -@kwdef struct Romberg <: IntegrationMethod - nsamples::Int = 1025 - stepsize::Float64 = 1.0 / 1024 - factors::Primes.Factorization{Int} = Primes.factor(1024) - nfactors::Int = 10 - atol::Float64 = 1e-6 +@kwdef struct Romberg{TI,TF,TP,TV} <: IntegrationMethod + nsamples::TI = 1025 + stepsize::TF = 1.0 / 1024 + factors::TP = Primes.factor(1024) + nfactors::TI = 10 + atol::TF = 1e-6 + samplepoints::TV = collect(range(0, 1, 1025)) end function set_romberg_options(; nsamples=1025, atol=1e-6) @@ -18,8 +19,8 @@ function set_romberg_options(; nsamples=1025, atol=1e-6) ) end -@kwdef struct GaussKronrod <: IntegrationMethod - order::Int = 3 - maxevals::Int = 1e2 - atol::Float64 = 1e-6 +@kwdef struct GaussKronrod{TI,TF} <: IntegrationMethod + order::TI = 3 + maxevals::TI = 1e2 + atol::TF = 1e-6 end From 0c4baee4e8d93fce09d3086f7ea255e31d41767b Mon Sep 17 00:00:00 2001 From: juddmehr Date: Thu, 28 Mar 2024 20:28:24 -0600 Subject: [PATCH 03/12] start implementing Gauss Legendre quadrature --- Project.toml | 1 + src/DuctAPE.jl | 5 +- src/precomputation/gausskronrod_integrals.jl | 5 +- src/precomputation/gausslegendre_integrals.jl | 77 +++++++++++++++++++ src/precomputation/romberg_integrals.jl | 7 +- src/types.jl | 12 ++- 6 files changed, 98 insertions(+), 9 deletions(-) create mode 100644 src/precomputation/gausslegendre_integrals.jl diff --git a/Project.toml b/Project.toml index 674fa3ab..43b68a8f 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.1.0" [deps] CCBlade = "e1828068-15df-11e9-03e4-ef195ea46fa4" FLOWMath = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3" +FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" ImplicitAD = "e7cbb90b-9b31-4eb2-a8c8-45099c074ee1" diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index 4a9bebbc..15919b72 100644 --- a/src/DuctAPE.jl +++ b/src/DuctAPE.jl @@ -11,7 +11,8 @@ include("C4Blade/C4Blade.jl") # augmented CCBlade implementation (cascade compat const c4b = C4Blade using SpecialFunctions # required for elliptic integrals -using QuadGK # required for integration of linear panels +using QuadGK # integration +using FastGaussQuadrature # integration using StaticArrays # used in miscellaneous places for code efficiency using LinearAlgebra: mul!, ldiv!, lu!, NoPivot, issuccess#, factorize # used in linear system assembly and solve @@ -54,6 +55,8 @@ include("precomputation/body_aic.jl") include("precomputation/induced_velocity_matrices.jl") include("precomputation/romberg_integrals.jl") +include("precomputation/gausslegendre_integrals.jl") +include("precomputation/gausskronrod_integrals.jl") include("precomputation/integrands.jl") ##### ----- SOLVER ----- ##### diff --git a/src/precomputation/gausskronrod_integrals.jl b/src/precomputation/gausskronrod_integrals.jl index 52da4960..c3c25fb8 100644 --- a/src/precomputation/gausskronrod_integrals.jl +++ b/src/precomputation/gausskronrod_integrals.jl @@ -18,14 +18,13 @@ function nominal_vortex_panel_integration( # Define function to integrate function fsample(t) - return nominal_vortex_induced_velocity_sample( + return dt.nominal_vortex_induced_velocity_sample( t, node1, node2, influence_length, controlpoint, - cache_vec; - nondimrange=nondimrange, + cache_vec ) end diff --git a/src/precomputation/gausslegendre_integrals.jl b/src/precomputation/gausslegendre_integrals.jl new file mode 100644 index 00000000..76cfcea0 --- /dev/null +++ b/src/precomputation/gausslegendre_integrals.jl @@ -0,0 +1,77 @@ +#---------------------------------# +# VORTEX # +#---------------------------------# + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function nominal_vortex_panel_integration!( + integration_options::GaussLegendre, + V, + node1, + node2, + influence_length, + controlpoint, + containers; + debug=false, +) + + # - Sample Function - # + for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) + # Define function to integrate + dt.nominal_vortex_induced_velocity_sample!( + s, t, node1, node2, influence_length, controlpoint, containers.sample_cache + ) + end + + # - Do Quadrature - # + for (i, s) in enumerate(eachcol(containers.samples)) + V[i] = dot(integration_options.weights, s) + end + + if debug + return V, err + else + return V + end +end + +#---------------------------------# +# SOURCE # +#---------------------------------# + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function nominal_source_panel_integration!( + integration_options::GaussLegendre, + V, + node1, + node2, + influence_length, + controlpoint, + containers; + debug=false, +) + + # - Sample Function - # + for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) + # Define function to integrate + # TODO: write this function + dt.nominal_source_induced_velocity_sample!( + s, t, node1, node2, influence_length, controlpoint, containers.sample_cache + ) + end + + # - Do Quadrature - # + for (i, s) in enumerate(eachcol(containers.samples)) + V[i] = dot(integration_options.weights, s) + end + + if debug + return V, err + else + return V + end +end + diff --git a/src/precomputation/romberg_integrals.jl b/src/precomputation/romberg_integrals.jl index 214ea9d0..5cfbcd3f 100644 --- a/src/precomputation/romberg_integrals.jl +++ b/src/precomputation/romberg_integrals.jl @@ -114,7 +114,6 @@ end #---------------------------------# # VORTEX # #---------------------------------# - """ `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ @@ -130,14 +129,14 @@ function nominal_vortex_panel_integration!( ) # Define function to integrate - for (s, t) in zip(eachrow(containers.samples), integration_options.samplepoints) - nominal_vortex_induced_velocity_sample!( + for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) + dt.nominal_vortex_induced_velocity_sample!( s, t, node1, node2, influence_length, controlpoint, containers.sample_cache ) end for (i, s) in enumerate(eachcol(containers.samples)) - V[i], containers.err[i] = romberg!( + V[i], containers.err[i] = dt.romberg!( containers.fh, integration_options.stepsize, s, diff --git a/src/types.jl b/src/types.jl index ba59725f..3c23041d 100644 --- a/src/types.jl +++ b/src/types.jl @@ -6,7 +6,7 @@ abstract type IntegrationMethod end factors::TP = Primes.factor(1024) nfactors::TI = 10 atol::TF = 1e-6 - samplepoints::TV = collect(range(0, 1, 1025)) + sample_points::TV = collect(range(0, 1, 1025)) end function set_romberg_options(; nsamples=1025, atol=1e-6) @@ -24,3 +24,13 @@ end maxevals::TI = 1e2 atol::TF = 1e-6 end + +struct GaussLegendre{TN,TW} <: IntegrationMethod + sample_points::TN + weights::TW +end + +function GaussLegendre(nsamples=20; silence_warnings=true) + nodes, weights = FastGaussQuadrature.gausslegendre(nsamples) + return GaussLegendre(linear_transform((-1, 1), (0, 1), nodes), weights ./ 2.0) +end From d1e34ceccf3416b1e3468a3d5efa594009c260fb Mon Sep 17 00:00:00 2001 From: juddmehr Date: Fri, 29 Mar 2024 21:20:11 -0600 Subject: [PATCH 04/12] get gauss legendre quadrature working for all cases --- src/precomputation/gausskronrod_integrals.jl | 2 +- src/precomputation/gausslegendre_integrals.jl | 87 +++- src/precomputation/integrals.jl | 370 +++++++++--------- src/precomputation/integrands.jl | 165 +++++++- src/types.jl | 12 +- 5 files changed, 435 insertions(+), 201 deletions(-) diff --git a/src/precomputation/gausskronrod_integrals.jl b/src/precomputation/gausskronrod_integrals.jl index c3c25fb8..396effed 100644 --- a/src/precomputation/gausskronrod_integrals.jl +++ b/src/precomputation/gausskronrod_integrals.jl @@ -18,7 +18,7 @@ function nominal_vortex_panel_integration( # Define function to integrate function fsample(t) - return dt.nominal_vortex_induced_velocity_sample( + return nominal_vortex_induced_velocity_sample( t, node1, node2, diff --git a/src/precomputation/gausslegendre_integrals.jl b/src/precomputation/gausslegendre_integrals.jl index 76cfcea0..546e4141 100644 --- a/src/precomputation/gausslegendre_integrals.jl +++ b/src/precomputation/gausslegendre_integrals.jl @@ -15,11 +15,10 @@ function nominal_vortex_panel_integration!( containers; debug=false, ) - # - Sample Function - # for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) # Define function to integrate - dt.nominal_vortex_induced_velocity_sample!( + nominal_vortex_induced_velocity_sample!( s, t, node1, node2, influence_length, controlpoint, containers.sample_cache ) end @@ -36,6 +35,47 @@ function nominal_vortex_panel_integration!( end end +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function self_vortex_panel_integration!( + integration_options::GaussLegendre, + V, + node1, + node2, + influence_length, + controlpoint, + containers; + debug=false, +) + + # - Sample Function - # + for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) + # Define function to integrate + self_vortex_induced_velocity_sample!( + s, t, node1, node2, influence_length, controlpoint, containers.sample_cache + ) + end + + # - Do Quadrature - # + for (i, s) in enumerate(eachcol(containers.samples)) + V[i] = dot(integration_options.weights, s) + end + + analytically_integrated_vortex_influence!( + @view(containers.sample_cache[1:2]), controlpoint[2], influence_length + ) + + V .*= influence_length + V[1:2] .+= containers.sample_cache[1] / 2.0 + + if debug + return reshape(V, (2, 2)), err + else + return reshape(V, (2, 2)) + end +end + #---------------------------------# # SOURCE # #---------------------------------# @@ -57,8 +97,7 @@ function nominal_source_panel_integration!( # - Sample Function - # for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) # Define function to integrate - # TODO: write this function - dt.nominal_source_induced_velocity_sample!( + nominal_source_induced_velocity_sample!( s, t, node1, node2, influence_length, controlpoint, containers.sample_cache ) end @@ -75,3 +114,43 @@ function nominal_source_panel_integration!( end end +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function self_source_panel_integration!( + integration_options::GaussLegendre, + V, + node1, + node2, + influence_length, + controlpoint, + containers; + debug=false, +) + + # - Sample Function - # + for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) + # Define function to integrate + self_source_induced_velocity_sample!( + s, t, node1, node2, influence_length, controlpoint, containers.sample_cache + ) + end + + # - Do Quadrature - # + for (i, s) in enumerate(eachcol(containers.samples)) + V[i] = dot(integration_options.weights, s) + end + + analytically_integrated_source_influence!( + @view(containers.sample_cache[1:2]), controlpoint[2], influence_length + ) + + V .*= influence_length + V[3:4] .+= containers.sample_cache[2] / 2.0 + + if debug + return reshape(V, (2, 2)), err + else + return reshape(V, (2, 2)) + end +end diff --git a/src/precomputation/integrals.jl b/src/precomputation/integrals.jl index 080b0204..64bb010b 100644 --- a/src/precomputation/integrals.jl +++ b/src/precomputation/integrals.jl @@ -1,185 +1,185 @@ -###################################################################### -# # -# Linear Vortex Panel Integration # -# # -###################################################################### - -#---------------------------------# -# Nominal # -#---------------------------------# - -""" -`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -""" -function nominal_vortex_panel_integration( - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=[0.0; 1.0], - debug=false, -) - - # Define function to integrate - function fsample(t) - return nominal_vortex_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, - ) - end - - V, err = quadgk(fsample, 0.0, 1.0; order=3, maxevals=1e2, atol=1e-6) - - if debug - return reshape(V, (2, 2)), err - # return V, err - else - return reshape(V, (2, 2)) - # return V - end -end - -#---------------------------------# -# Singular # -#---------------------------------# - -""" -`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -""" -function self_vortex_panel_integration( - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=[0.0; 1.0], - debug=false, -) - - # Define function to integrate - function fsample(t) - return self_vortex_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, - ) - end - - V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) - - cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( - controlpoint[2], influence_length - ) - - V .*= influence_length - V[1:2] .+= cache_vec[1] / 2.0 - - if debug - return reshape(V, (2, 2)), err - else - return reshape(V, (2, 2)) - end -end - -###################################################################### -# # -# Linear Source Panel Integration # -# # -###################################################################### - -#---------------------------------# -# Nominal # -#---------------------------------# - -""" -`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -""" -function nominal_source_panel_integration( - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=[0.0; 1.0], - debug=false, -) - - # Define function to integrate - function fsample(t) - return nominal_source_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, - ) - end - - V, err = quadgk(fsample, 0.0, 1.0; order=3, maxevals=1e2, atol=1e-6) - - if debug - return reshape(V, (2, 2)), err - # return V, err - else - return reshape(V, (2, 2)) - # return V - end -end - -#---------------------------------# -# Singular # -#---------------------------------# - -""" -`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -""" -function self_source_panel_integration( - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=[0.0; 1.0], - debug=false, -) - - # Define function to integrate - function fsample(t) - return self_source_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, - ) - end - - V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) - - cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( - controlpoint[2], influence_length - ) - - V .*= influence_length - V[3:4] .+= cache_vec[2] / 2.0 - - if debug - return reshape(V, (2, 2)), err - else - return reshape(V, (2, 2)) - end -end +####################################################################### +## # +## Linear Vortex Panel Integration # +## # +####################################################################### + +##---------------------------------# +## Nominal # +##---------------------------------# + +#""" +#`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +#""" +#function nominal_vortex_panel_integration( +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=[0.0; 1.0], +# debug=false, +#) + +# # Define function to integrate +# function fsample(t) +# return nominal_vortex_induced_velocity_sample( +# t, +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=nondimrange, +# ) +# end + +# V, err = quadgk(fsample, 0.0, 1.0; order=3, maxevals=1e2, atol=1e-6) + +# if debug +# return reshape(V, (2, 2)), err +# # return V, err +# else +# return reshape(V, (2, 2)) +# # return V +# end +#end + +##---------------------------------# +## Singular # +##---------------------------------# + +#""" +#`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +#""" +#function self_vortex_panel_integration( +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=[0.0; 1.0], +# debug=false, +#) + +# # Define function to integrate +# function fsample(t) +# return self_vortex_induced_velocity_sample( +# t, +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=nondimrange, +# ) +# end + +# V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) + +# cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( +# controlpoint[2], influence_length +# ) + +# V .*= influence_length +# V[1:2] .+= cache_vec[1] / 2.0 + +# if debug +# return reshape(V, (2, 2)), err +# else +# return reshape(V, (2, 2)) +# end +#end + +####################################################################### +## # +## Linear Source Panel Integration # +## # +####################################################################### + +##---------------------------------# +## Nominal # +##---------------------------------# + +#""" +#`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +#""" +#function nominal_source_panel_integration( +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=[0.0; 1.0], +# debug=false, +#) + +# # Define function to integrate +# function fsample(t) +# return nominal_source_induced_velocity_sample( +# t, +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=nondimrange, +# ) +# end + +# V, err = quadgk(fsample, 0.0, 1.0; order=3, maxevals=1e2, atol=1e-6) + +# if debug +# return reshape(V, (2, 2)), err +# # return V, err +# else +# return reshape(V, (2, 2)) +# # return V +# end +#end + +##---------------------------------# +## Singular # +##---------------------------------# + +#""" +#`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +#""" +#function self_source_panel_integration( +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=[0.0; 1.0], +# debug=false, +#) + +# # Define function to integrate +# function fsample(t) +# return self_source_induced_velocity_sample( +# t, +# node1, +# node2, +# influence_length, +# controlpoint, +# cache_vec; +# nondimrange=nondimrange, +# ) +# end + +# V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) + +# cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( +# controlpoint[2], influence_length +# ) + +# V .*= influence_length +# V[3:4] .+= cache_vec[2] / 2.0 + +# if debug +# return reshape(V, (2, 2)), err +# else +# return reshape(V, (2, 2)) +# end +#end diff --git a/src/precomputation/integrands.jl b/src/precomputation/integrands.jl index 65275aff..449bc8d9 100644 --- a/src/precomputation/integrands.jl +++ b/src/precomputation/integrands.jl @@ -17,8 +17,8 @@ function nominal_vortex_induced_velocity_sample!( ) # Transform from (0,1) to actual position on panel - cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) # z coordinate - cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) # r coordinate + cache_vec[1] = linear_transform(nondimrange, (node1[1], node2[1]), t) # z coordinate + cache_vec[2] = linear_transform(nondimrange, (node1[2], node2[2]), t) # r coordinate # get relative geometry: xi, rho, m, r0 = calculate_xrm(controlpoint, [z; r]) calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) @@ -60,8 +60,8 @@ function nominal_vortex_induced_velocity_sample( ) # Transform from (0,1) to actual position on panel - cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) # z coordinate - cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) # r coordinate + cache_vec[1] = linear_transform(nondimrange, (node1[1], node2[1]), t) # z coordinate + cache_vec[2] = linear_transform(nondimrange, (node1[2], node2[2]), t) # r coordinate # get relative geometry: xi, rho, m, r0 = calculate_xrm(controlpoint, [z; r]) calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) @@ -152,6 +152,14 @@ function subtracted_singular_vortex_influence!(node, controlpoint, cache_vec) return axial, radial end +""" +""" +function analytically_integrated_vortex_influence!(V, r, influence_length) + V[1] = (influence_length / (4.0 * pi * r)) * (1.0 + log(16.0 * r / influence_length)) + V[2] = 0.0 + return V +end + """ """ function analytically_integrated_vortex_influence(r, influence_length) @@ -161,7 +169,52 @@ function analytically_integrated_vortex_influence(r, influence_length) end """ +# Arguments: +- `t::Float` : sample point in range (0,1) selected by quadrature. +""" +function self_vortex_induced_velocity_sample!( + V, t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] +) + + # Transform from (0,1) to actual position on panel + cache_vec[1] = linear_transform(nondimrange, (node1[1], node2[1]), t) + cache_vec[2] = linear_transform(nondimrange, (node1[2], node2[2]), t) + # get relative geometry + calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) + + # Get velocity components at sample points + cache_vec[7] = vortex_ring_vz!( + cache_vec[3],#xi + cache_vec[4],#rho + cache_vec[5],#m + cache_vec[2],#rj + 1.0, + view(cache_vec, 11:15), + ) #vz + cache_vec[8] = vortex_ring_vr!( + cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[2], view(cache_vec, 11:16) + ) #vr + + # Get singular piece to subtract + cache_vec[9], cache_vec[10] = subtracted_singular_vortex_influence!( + (cache_vec[1], cache_vec[2]), controlpoint, view(cache_vec, 11:15) + ) + + #= + assemble output components in the format: + [x_j r_j; x_{j+1} r_{j+1}] + and scale by influence panel length + (due to transformation of integration range to/from range=(0,1)) + =# + V[1] = cache_vec[7] * (1.0 - t) - cache_vec[9] / 2.0 + V[2] = cache_vec[7] * t .- cache_vec[9] / 2.0 + V[3] = cache_vec[8] * (1.0 - t) .- cache_vec[10] / 2.0 + V[4] = cache_vec[8] * t .- cache_vec[10] / 2.0 + + return V +end +""" # Arguments: - `t::Float` : sample point in range (0,1) selected by quadrature. """ @@ -170,8 +223,8 @@ function self_vortex_induced_velocity_sample( ) # Transform from (0,1) to actual position on panel - cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) - cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) + cache_vec[1] = linear_transform(nondimrange, (node1[1], node2[1]), t) + cache_vec[2] = linear_transform(nondimrange, (node1[2], node2[2]), t) # get relative geometry calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) @@ -219,7 +272,46 @@ end #---------------------------------# """ +# Arguments: +- `t::Float` : sample point in range (0,1) selected by quadrature. +""" +function nominal_source_induced_velocity_sample!( + V, t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] +) + + # Transform from (0,1) to actual position on panel + cache_vec[1] = linear_transform(nondimrange, (node1[1], node2[1]), t) + cache_vec[2] = linear_transform(nondimrange, (node1[2], node2[2]), t) + + # get relative geometry + calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) + + # Get velocity components at sample points + cache_vec[7] = source_ring_vz!( + cache_vec[3],#xi + cache_vec[4],#rho + cache_vec[5],#m + cache_vec[6],#rj + view(cache_vec, 11:15), + ) + cache_vec[8] = source_ring_vr!( + cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[6], view(cache_vec, 11:16) + ) + + #= + assemble output components in the format: + [x_j r_j; x_{j+1} r_{j+1}] + and scale by influence panel length + (due to transformation of integration range to/from range=(0,1)) + =# + V[1] = cache_vec[7] * (1.0 - t) * influence_length + V[2] = cache_vec[7] * t * influence_length + V[3] = cache_vec[8] * (1.0 - t) * influence_length + V[4] = cache_vec[8] * t * influence_length + return V +end +""" # Arguments: - `t::Float` : sample point in range (0,1) selected by quadrature. """ @@ -228,8 +320,8 @@ function nominal_source_induced_velocity_sample( ) # Transform from (0,1) to actual position on panel - cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) - cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) + cache_vec[1] = linear_transform(nondimrange, (node1[1], node2[1]), t) + cache_vec[2] = linear_transform(nondimrange, (node1[2], node2[2]), t) # get relative geometry calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) @@ -301,6 +393,14 @@ function subtracted_singular_source_influence!(node, controlpoint, cache_vec) return axial, radial end +""" +""" +function analytically_integrated_source_influence!(V, r, influence_length) + V[2] = (influence_length / (4.0 * pi * r)) * (1.0 + log(2.0 * r / influence_length)) + V[1] = 0.0 + return V +end + """ """ function analytically_integrated_source_influence(r, influence_length) @@ -310,7 +410,52 @@ function analytically_integrated_source_influence(r, influence_length) end """ +# Arguments: +- `t::Float` : sample point in range (0,1) selected by quadrature. +""" +function self_source_induced_velocity_sample!( + V, t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] +) + # Transform from (0,1) to actual position on panel + cache_vec[1] = linear_transform(nondimrange, (node1[1], node2[1]), t) + cache_vec[2] = linear_transform(nondimrange, (node1[2], node2[2]), t) + + # get relative geometry + calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) + + # Get velocity components at sample points + cache_vec[7] = source_ring_vz!( + cache_vec[3],#xi + cache_vec[4],#rho + cache_vec[5],#m + cache_vec[2],#rj + view(cache_vec, 11:15), + ) + cache_vec[8] = source_ring_vr!( + cache_vec[3], cache_vec[4], cache_vec[5], cache_vec[2], view(cache_vec, 11:16) + ) + + # Get singular piece to subtract + cache_vec[9], cache_vec[10] = subtracted_singular_source_influence!( + (cache_vec[1], cache_vec[2]), controlpoint, view(cache_vec, 11:15) + ) + + #= + assemble output components in the format: + [x_j r_j; x_{j+1} r_{j+1}] + and scale by influence panel length + (due to transformation of integration range to/from range=(0,1)) + =# + V[1] = cache_vec[7] * (1.0 - t) .- cache_vec[9] / 2.0 + V[2] = cache_vec[7] * t .- cache_vec[9] / 2.0 + V[3] = cache_vec[8] * (1.0 - t) .- cache_vec[10] / 2.0 + V[4] = cache_vec[8] * t .- cache_vec[10] / 2.0 + + return V +end + +""" # Arguments: - `t::Float` : sample point in range (0,1) selected by quadrature. """ @@ -319,8 +464,8 @@ function self_source_induced_velocity_sample( ) # Transform from (0,1) to actual position on panel - cache_vec[1] = fm.linear(nondimrange, (node1[1], node2[1]), t) - cache_vec[2] = fm.linear(nondimrange, (node1[2], node2[2]), t) + cache_vec[1] = linear_transform(nondimrange, (node1[1], node2[1]), t) + cache_vec[2] = linear_transform(nondimrange, (node1[2], node2[2]), t) # get relative geometry calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) diff --git a/src/types.jl b/src/types.jl index 3c23041d..ba010348 100644 --- a/src/types.jl +++ b/src/types.jl @@ -21,7 +21,7 @@ end @kwdef struct GaussKronrod{TI,TF} <: IntegrationMethod order::TI = 3 - maxevals::TI = 1e2 + maxevals::TI = 100 atol::TF = 1e-6 end @@ -31,6 +31,16 @@ struct GaussLegendre{TN,TW} <: IntegrationMethod end function GaussLegendre(nsamples=20; silence_warnings=true) + if silence_warnings && Bool((nsamples) % 2) + @warn "Must have an even number of GaussLegendre sample points if using for panel self influence" + end + nodes, weights = FastGaussQuadrature.gausslegendre(nsamples) + return GaussLegendre(linear_transform((-1, 1), (0, 1), nodes), weights ./ 2.0) end + +@kwdef struct IntegrationOptions{TN,TS} + nominal::TN = GaussLegendre(20) + singular::TS = GaussKronrod(3, 1e2, 1e-6) +end From f36dcf82904ac982642da53845ceb23d69634486 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Mon, 1 Apr 2024 07:52:29 -0600 Subject: [PATCH 05/12] get integration drafted for all cases for romberg integration. --- src/precomputation/gausskronrod_integrals.jl | 38 +- src/precomputation/integrands.jl | 32 +- src/precomputation/romberg_integrals.jl | 430 +++++++++++-------- src/types.jl | 24 +- 4 files changed, 273 insertions(+), 251 deletions(-) diff --git a/src/precomputation/gausskronrod_integrals.jl b/src/precomputation/gausskronrod_integrals.jl index 396effed..e76f3090 100644 --- a/src/precomputation/gausskronrod_integrals.jl +++ b/src/precomputation/gausskronrod_integrals.jl @@ -12,19 +12,13 @@ function nominal_vortex_panel_integration( influence_length, controlpoint, cache_vec; - nondimrange=[0.0; 1.0], debug=false, ) # Define function to integrate function fsample(t) return nominal_vortex_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec + t, node1, node2, influence_length, controlpoint, cache_vec ) end @@ -50,25 +44,19 @@ end `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ function self_vortex_panel_integration( + integration_options::GaussKronrod, node1, node2, influence_length, controlpoint, cache_vec; - nondimrange=[0.0; 1.0], debug=false, ) # Define function to integrate function fsample(t) return self_vortex_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, + t, node1, node2, influence_length, controlpoint, cache_vec; ) end @@ -104,25 +92,19 @@ end `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ function nominal_source_panel_integration( + integration_options::GaussKronrod, node1, node2, influence_length, controlpoint, cache_vec; - nondimrange=[0.0; 1.0], debug=false, ) # Define function to integrate function fsample(t) return nominal_source_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, + t, node1, node2, influence_length, controlpoint, cache_vec; ) end @@ -148,25 +130,19 @@ end `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ function self_source_panel_integration( + integration_options::GaussKronrod, node1, node2, influence_length, controlpoint, cache_vec; - nondimrange=[0.0; 1.0], debug=false, ) # Define function to integrate function fsample(t) return self_source_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, + t, node1, node2, influence_length, controlpoint, cache_vec; ) end diff --git a/src/precomputation/integrands.jl b/src/precomputation/integrands.jl index 449bc8d9..4511a836 100644 --- a/src/precomputation/integrands.jl +++ b/src/precomputation/integrands.jl @@ -43,10 +43,10 @@ function nominal_vortex_induced_velocity_sample!( and scale by influence panel length (due to transformation of integration range to/from range=(0,1)) =# - V[1] = cache_vec[7] * (1.0 - t) * influence_length - V[2] = cache_vec[7] * t * influence_length - V[3] = cache_vec[8] * (1.0 - t) * influence_length - V[4] = cache_vec[8] * t * influence_length + V[1] += cache_vec[7] * (1.0 - t) * influence_length + V[2] += cache_vec[7] * t * influence_length + V[3] += cache_vec[8] * (1.0 - t) * influence_length + V[4] += cache_vec[8] * t * influence_length return V end @@ -207,10 +207,10 @@ function self_vortex_induced_velocity_sample!( and scale by influence panel length (due to transformation of integration range to/from range=(0,1)) =# - V[1] = cache_vec[7] * (1.0 - t) - cache_vec[9] / 2.0 - V[2] = cache_vec[7] * t .- cache_vec[9] / 2.0 - V[3] = cache_vec[8] * (1.0 - t) .- cache_vec[10] / 2.0 - V[4] = cache_vec[8] * t .- cache_vec[10] / 2.0 + V[1] += cache_vec[7] * (1.0 - t) - cache_vec[9] / 2.0 + V[2] += cache_vec[7] * t .- cache_vec[9] / 2.0 + V[3] += cache_vec[8] * (1.0 - t) .- cache_vec[10] / 2.0 + V[4] += cache_vec[8] * t .- cache_vec[10] / 2.0 return V end @@ -304,10 +304,10 @@ function nominal_source_induced_velocity_sample!( and scale by influence panel length (due to transformation of integration range to/from range=(0,1)) =# - V[1] = cache_vec[7] * (1.0 - t) * influence_length - V[2] = cache_vec[7] * t * influence_length - V[3] = cache_vec[8] * (1.0 - t) * influence_length - V[4] = cache_vec[8] * t * influence_length + V[1] += cache_vec[7] * (1.0 - t) * influence_length + V[2] += cache_vec[7] * t * influence_length + V[3] += cache_vec[8] * (1.0 - t) * influence_length + V[4] += cache_vec[8] * t * influence_length return V end @@ -447,10 +447,10 @@ function self_source_induced_velocity_sample!( and scale by influence panel length (due to transformation of integration range to/from range=(0,1)) =# - V[1] = cache_vec[7] * (1.0 - t) .- cache_vec[9] / 2.0 - V[2] = cache_vec[7] * t .- cache_vec[9] / 2.0 - V[3] = cache_vec[8] * (1.0 - t) .- cache_vec[10] / 2.0 - V[4] = cache_vec[8] * t .- cache_vec[10] / 2.0 + V[1] += cache_vec[7] * (1.0 - t) .- cache_vec[9] / 2.0 + V[2] += cache_vec[7] * t .- cache_vec[9] / 2.0 + V[3] += cache_vec[8] * (1.0 - t) .- cache_vec[10] / 2.0 + V[4] += cache_vec[8] * t .- cache_vec[10] / 2.0 return V end diff --git a/src/precomputation/romberg_integrals.jl b/src/precomputation/romberg_integrals.jl index 5cfbcd3f..f8dfd4da 100644 --- a/src/precomputation/romberg_integrals.jl +++ b/src/precomputation/romberg_integrals.jl @@ -1,12 +1,3 @@ -#TODO: decide if dfdc method or julia method is better -# the DFDC method requires the whole extrapolation, but you can limit the number of intervals to evaluate when the integral is easy. -# when the integral is hard though, you end up doing 2x more calculations than you would have if you set your sample size larger to begin with. -# the julia method requires you to evaluate all the sample points up front, but will terminate the extrapolation as soon as you hit your desired tolerance. So you need to have a large number of samples to start with. -# DFDC doesn't actually do a true trapezoidal integration, rather it assumes a linear function and does δx*f(midpoint of iterval), this is technically less accurate than the trapezoidal method, but does avoid sampling at the end and midpoints of panels. Theoretically, the extrapolation will account for the less accurate integration method. -# perhaps a combo, where you do a DFDC-like integration (not quite trapezoidal), but do the large sample size and early extrapolation termination. - -# Current Plan: do romberg for the non-singular ones, and quadgk for the singular ones. - """ extrapolate!(fh::AbstractVector; power=1, atol=0, rtol=0, maxeval=typemax(Int), breaktol=Inf) @@ -19,104 +10,76 @@ the type of `h` and `T` is the type of the extrapolated `f(0)` **result**. This should be a floating-point type, i.e. `fh` should contain `float(f(h))` if the function you are extrapolating is not already floating-point-valued.) """ -function extrapolate!( - fh::AbstractVector{<:Tuple{TF1,TF2}}; - power::Int=2, - atol::Float64=1e-6, - rtol::Float64=0.0, - breaktol::Float64=Inf, - maxeval::Int=typemax(Int), -) where {TF1,TF2} - - # - Check for Problems - # - (rtol >= 0 && atol >= zero(atol)) || - throw(ArgumentError("rtol and atol must be nonnegative")) - breaktol > 0 || throw(ArgumentError("breaktol must be positive")) - isempty(fh) && throw(ArgumentError("(f,h) array must be non-empty")) +function extrapolate!(V, err, fh; power=2, atol=1e-6) where {TF1,TF2} # - rename for convenience - # (f0, h) = first(fh) # - initialize - # - err::typeof(float(norm(f0))) = Inf + err .= Inf numeval = 1 - maxeval = min(maxeval, length(fh)) + maxeval = size(fh, 1) # - Loop - # while numeval < maxeval + # - increment - # numeval += 1 # - update - # (f_prime, h_prime) = fh[numeval + (firstindex(fh) - 1)] - abs(h) > abs(h_prime) || - throw(ArgumentError("|$h_prime| >= |$h| is not decreasing")) h = h_prime - minerr_prime = oftype(err, Inf) + minerr_prime = similar(err) .= Inf f_ip1 = f_prime # - Inner Loop - # for i in (numeval - 1):-1:1 f_i, h_i = fh[i + (firstindex(fh) - 1)] - c = (h_i / h_prime)^power - f_ip1 += (f_ip1 - f_i) / (c - 1) + c = (h_i[] / h_prime[])^power + println((f_ip1 - f_i)) + @. f_ip1 += (f_ip1 - f_i) / (c - 1) fh[i + (firstindex(fh) - 1)] = (f_ip1, h_i) - err_prime = norm(f_ip1 - f_i) - minerr_prime = min(minerr_prime, err_prime) + err_prime = norm.(f_ip1 - f_i) + minerr_prime = min.(minerr_prime, err_prime) if err_prime < err - f0, err = f_ip1, err_prime + V, err = f_ip1, err_prime end end - # - stop early if error increases too much - # - (minerr_prime > breaktol * err || !isfinite(minerr_prime)) && break - # - converged - # - err <= max(rtol * norm(f0), atol) && break + all(err .<= atol) && break end - return (f0, err) -end -@views function romberg!(fh, step_size, samples, endsum, factors, numfactors; kwargs...) - - # - Do Trapezoidal Integration with Composite Trapezoidal Rule - # - b, e = firstindex(samples), lastindex(samples) - i = numfactors + 1 - fh[i] = (step_size * (sum(samples[(b + 1):(e - 1)]) + endsum), step_size) - nstep = 1 - for (j, K) in factors - @inbounds for k in 1:K - nstep *= j - sdx = nstep * step_size - if i == 2 # last iteration (empty sum) - fh[1] = (sdx * endsum, sdx) - else - fh[i -= 1] = ( - sdx * (sum(samples[(b + nstep):nstep:(e - nstep)]) + endsum), sdx - ) - end - end - end + return V, err +end - # - Do Richardson Extrapolation until converged - # - return extrapolate!(fh; power=2, kwargs...) +function allocate_integration_containers( + integration_options::Romberg, dispatch_type::TF; cache_size=20 +) where {TF} + return (; + sample_cache=zeros(TF, cache_size), + samples=[(zeros(TF, 4), TF[0.0]) for i in 1:(integration_options.max_subdivisions)], + ) end -function allocate_integration_containers(dispatch_type::TF, nsamples, nfactors) where {TF} +function allocate_integration_containers( + integration_options::GaussLegendre, dispatch_type::TF; cache_size=20 +) where {TF} return (; - err=zeros(TF, 2, 2), #always this size - sample_cache=zeros(TF, 20), #always this size, this is for the intermediate calcs - samples=zeros(TF, nsamples, 4), # always 4 wide - fh=[(TF(1.0), TF(1.0)) for i in 1:(nfactors + 1)], + sample_cache=zeros(TF, cache_size), samples=zeros(TF, integration_options.nsamples) ) end +function allocate_integration_containers( + integration_options::GaussKronrod, dispatch_type::TF; cache_size=20 +) where {TF} + return (; sample_cache=zeros(TF, cache_size)) +end + #---------------------------------# # VORTEX # #---------------------------------# -""" -`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -""" + function nominal_vortex_panel_integration!( integration_options::Romberg, V, @@ -128,155 +91,248 @@ function nominal_vortex_panel_integration!( debug=false, ) - # Define function to integrate - for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) - dt.nominal_vortex_induced_velocity_sample!( - s, t, node1, node2, influence_length, controlpoint, containers.sample_cache - ) + # - Loop through number of subdivisions (start with 2) - # + for ii in 1:(integration_options.max_subdivisions) + nint = 2^ii + + # - Set step size for integration and extrapolation - # + dx = 1.0 / nint + containers.samples[ii][2] .= dx + + # - Loop through intervals for current subdivision level - # + for i in 1:nint + # sample at the interval midpoints + t = (i - 0.5) / nint + + # get sample velocity + dt.nominal_vortex_induced_velocity_sample!( + @view(containers.samples[ii])[1][1], + t, + node1, + node2, + influence_length, + controlpoint, + containers.sample_cache, + ) + end + + # multiply samples by the interval length + containers.samples[ii][1] .*= dx + + if ii > 1 + # extrapolate once you can + dt.extrapolate!( + V, + @view(containers.sample_cache[1:4]), + @view(containers.samples[1:ii]); + atol=integration_options.atol, + ) + + # if you have met the convergence requirements, terminate + all(@view(containers.sample_cache[1:4]) .<= integration_options.atol) && break + end end - for (i, s) in enumerate(eachcol(containers.samples)) - V[i], containers.err[i] = dt.romberg!( - containers.fh, - integration_options.stepsize, - s, - (s[1] + s[end]) / 2.0, - integration_options.factors, - integration_options.nfactors; - atol=integration_options.atol, + return V +end + +""" +`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] +""" +function self_vortex_panel_integration( + integration_options::Romberg, + node1, + node2, + influence_length, + controlpoint, + containers, + debug=false, +) + + # - Loop through number of subdivisions (start with 2) - # + for ii in 1:(integration_options.max_subdivisions) + nint = 2^ii + + # - Set step size for integration and extrapolation - # + dx = 1.0 / nint + containers.samples[ii][2] .= dx + + # - Loop through intervals for current subdivision level - # + for i in 1:nint + # sample at the interval midpoints + t = (i - 0.5) / nint + + # get sample velocity + dt.self_vortex_induced_velocity_sample!( + @view(containers.samples[ii])[1][1], + t, + node1, + node2, + influence_length, + controlpoint, + containers.sample_cache, + ) + end + + # multiply samples by the interval length + containers.samples[ii][1] .*= dx + # multiply samples by the influence length (not done in self induced unit velocity function) + containers.samples[ii][1] .*= influence_length + + # - Add in the analytical bits - # + containers.sample_cache[1], containers.sample_cache[2] = analytically_integrated_vortex_influence( + controlpoint[2], influence_length ) + containers.samples[ii][1][1:2] .+= containers.sample_cache[1] / 2.0 + + if ii > 1 + # extrapolate once you can + dt.extrapolate!( + V, + @view(containers.sample_cache[1:4]), + @view(containers.samples[1:ii]); + atol=integration_options.atol, + ) + + # if you have met the convergence requirements, terminate + all(@view(containers.sample_cache[1:4]) .<= integration_options.atol) && break + end end if debug - # return reshape(V, (2, 2)), err return V, err else - # return reshape(V, (2, 2)) return V end end -# """ -# `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -# """ -# function self_vortex_panel_integration( -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=[0.0; 1.0], -# debug=false, -# ) - -# # Define function to integrate -# function fsample(t) -# return self_vortex_induced_velocity_sample( -# t, -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=nondimrange, -# ) -# end - -# V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) - -# cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( -# controlpoint[2], influence_length -# ) - -# V .*= influence_length -# V[1:2] .+= cache_vec[1] / 2.0 - -# if debug -# return reshape(V, (2, 2)), err -# else -# return reshape(V, (2, 2)) -# end -# end - #---------------------------------# # SOURCE # #---------------------------------# +function nominal_source_panel_integration!( + integration_options::Romberg, + V, + node1, + node2, + influence_length, + controlpoint, + containers; + debug=false, +) + + # - Loop through number of subdivisions (start with 2) - # + for ii in 1:(integration_options.max_subdivisions) + nint = 2^ii + + # - Set step size for integration and extrapolation - # + dx = 1.0 / nint + containers.samples[ii][2] .= dx + + # - Loop through intervals for current subdivision level - # + for i in 1:nint + # sample at the interval midpoints + t = (i - 0.5) / nint + + # get sample velocity + dt.nominal_source_induced_velocity_sample!( + @view(containers.samples[ii])[1][1], + t, + node1, + node2, + influence_length, + controlpoint, + containers.sample_cache, + ) + end + + # multiply samples by the interval length + containers.samples[ii][1] .*= dx + + if ii > 1 + # extrapolate once you can + dt.extrapolate!( + V, + @view(containers.sample_cache[1:4]), + @view(containers.samples[1:ii]); + atol=integration_options.atol, + ) + + # if you have met the convergence requirements, terminate + all(@view(containers.sample_cache[1:4]) .<= integration_options.atol) && break + end + end + + return V +end + """ `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ -function nominal_source_panel_integration!( +function self_source_panel_integration( + integration_options::Romberg, node1, node2, influence_length, controlpoint, - cache_vec; - nondimrange=[0.0; 1.0], + containers; debug=false, ) - # Define function to integrate - function fsample(t) - return nominal_source_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, + # - Loop through number of subdivisions (start with 2) - # + for ii in 1:(integration_options.max_subdivisions) + nint = 2^ii + + # - Set step size for integration and extrapolation - # + dx = 1.0 / nint + containers.samples[ii][2] .= dx + + # - Loop through intervals for current subdivision level - # + for i in 1:nint + # sample at the interval midpoints + t = (i - 0.5) / nint + + # get sample velocity + dt.self_source_induced_velocity_sample!( + @view(containers.samples[ii])[1][1], + t, + node1, + node2, + influence_length, + controlpoint, + containers.sample_cache, + ) + end + + # multiply samples by the interval length + containers.samples[ii][1] .*= dx + # multiply samples by the influence length (not done in self induced unit velocity function) + containers.samples[ii][1] .*= influence_length + + # - Add in the analytical bits - # + containers.sample_cache[1], containers.sample_cache[2] = analytically_integrated_source_influence( + controlpoint[2], influence_length ) + containers.samples[ii][1][3:4] .+= containers.sample_cache[2] / 2.0 + + if ii > 1 + # extrapolate once you can + dt.extrapolate!( + V, + @view(containers.sample_cache[1:4]), + @view(containers.samples[1:ii]); + atol=integration_options.atol, + ) + + # if you have met the convergence requirements, terminate + all(@view(containers.sample_cache[1:4]) .<= integration_options.atol) && break + end end - V, err = quadgk(fsample, 0.0, 1.0; order=3, maxevals=1e2, atol=1e-6) - if debug - return reshape(V, (2, 2)), err - # return V, err + return V, err else - return reshape(V, (2, 2)) - # return V + return V end end -# """ -# `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -# """ -# function self_source_panel_integration( -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=[0.0; 1.0], -# debug=false, -# ) - -# # Define function to integrate -# function fsample(t) -# return self_source_induced_velocity_sample( -# t, -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=nondimrange, -# ) -# end - -# V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) - -# cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( -# controlpoint[2], influence_length -# ) - -# V .*= influence_length -# V[3:4] .+= cache_vec[2] / 2.0 - -# if debug -# return reshape(V, (2, 2)), err -# else -# return reshape(V, (2, 2)) -# end -# end diff --git a/src/types.jl b/src/types.jl index ba010348..325b95cd 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,28 +1,18 @@ abstract type IntegrationMethod end @kwdef struct Romberg{TI,TF,TP,TV} <: IntegrationMethod - nsamples::TI = 1025 - stepsize::TF = 1.0 / 1024 - factors::TP = Primes.factor(1024) - nfactors::TI = 10 + max_subdivisions::TI = 10 atol::TF = 1e-6 - sample_points::TV = collect(range(0, 1, 1025)) end -function set_romberg_options(; nsamples=1025, atol=1e-6) - stepsize = 1.0 / (nsamples - 1) - factors = Primes.factor(nsamples - 1) - nfactors = sum(values(factors)) - - return Romberg(; - atol=atol, nsamples=nsamples, factors=factors, nfactors=nfactors, stepsize=stepsize - ) +function set_romberg_options(; max_subdivisions=10, atol=1e-6) + return Romberg(; max_subdivisions, atol) end @kwdef struct GaussKronrod{TI,TF} <: IntegrationMethod - order::TI = 3 - maxevals::TI = 100 - atol::TF = 1e-6 + order::TI = 5 + maxevals::TI = 1000 + atol::TF = 1e-12 end struct GaussLegendre{TN,TW} <: IntegrationMethod @@ -42,5 +32,5 @@ end @kwdef struct IntegrationOptions{TN,TS} nominal::TN = GaussLegendre(20) - singular::TS = GaussKronrod(3, 1e2, 1e-6) + singular::TS = GaussLegendre(20) end From b3f64bf6696dcd85f10e90c9eb5b47d74f89b056 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Wed, 3 Apr 2024 19:23:20 -0600 Subject: [PATCH 06/12] fix typo --- src/preprocess/velocities/integrals.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/preprocess/velocities/integrals.jl b/src/preprocess/velocities/integrals.jl index 9c45e67d..64bb010b 100644 --- a/src/preprocess/velocities/integrals.jl +++ b/src/preprocess/velocities/integrals.jl @@ -1,4 +1,3 @@ -<<<<<<< HEAD:src/precomputation/integrals.jl ####################################################################### ## # ## Linear Vortex Panel Integration # From 7da2a4f397f94be6e8e6b2494c4b1c76129f5109 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Wed, 3 Apr 2024 19:30:40 -0600 Subject: [PATCH 07/12] move over types --- src/DuctAPE.jl | 18 +++++++++--------- src/types.jl | 36 ------------------------------------ src/utilities/types.jl | 40 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 49 insertions(+), 45 deletions(-) diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index a86fa466..3b62e3d0 100644 --- a/src/DuctAPE.jl +++ b/src/DuctAPE.jl @@ -12,6 +12,7 @@ const c4b = C4Blade # - Packages for Calculating Unit Induced Velocities - # # For Kernels using SpecialFunctions # required for elliptic integrals + # For Integration using FastGaussQuadrature using QuadGK @@ -23,11 +24,14 @@ using PreallocationTools # caches # - Packages for Solves - # using ImplicitAD # used for all solves using LinearAlgebra # linear solve and LU decomposition + # General Nonlinear solves using SimpleNonlinearSolve + # Fixed-Point Iteration Solvers # using FixedPointAcceleration using SpeedMapping + # For using NLsolve using NLsolve #for newton solver using LineSearches # used in newton solver @@ -88,24 +92,20 @@ include("preprocess/geometry/elliptic_grid_residuals.jl") # Induced Velocity Functions include("preprocess/velocities/unit_induced_velocities.jl") -include("preprocess/velocities/integrals.jl") include("preprocess/velocities/induced_velocity_matrices.jl") include("preprocess/velocities/body_aic.jl") +# Quadrature +include("preprocess/velocities/integrals.jl") +include("preprocess/velocities/gausslegendre_integrals.jl") +include("preprocess/velocities/gausskronrod_integrals.jl") +include("preprocess/velocities/integrands.jl") ##### ----- PROCESS ----- ##### - # Solve and Residual Functions include("process/solve.jl") include("process/residuals/CSORresidual.jl") include("process/residuals/systemresidual.jl") -include("precomputation/gausslegendre_integrals.jl") -include("precomputation/gausskronrod_integrals.jl") -include("precomputation/integrands.jl") -##### ----- SOLVER ----- ##### - -include("solve/solve.jl") - # Aerodynamics Functions include("process/aerodynamics/body_aerodynamics.jl") include("process/aerodynamics/rotor_aerodynamics.jl") diff --git a/src/types.jl b/src/types.jl index 325b95cd..e69de29b 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,36 +0,0 @@ -abstract type IntegrationMethod end - -@kwdef struct Romberg{TI,TF,TP,TV} <: IntegrationMethod - max_subdivisions::TI = 10 - atol::TF = 1e-6 -end - -function set_romberg_options(; max_subdivisions=10, atol=1e-6) - return Romberg(; max_subdivisions, atol) -end - -@kwdef struct GaussKronrod{TI,TF} <: IntegrationMethod - order::TI = 5 - maxevals::TI = 1000 - atol::TF = 1e-12 -end - -struct GaussLegendre{TN,TW} <: IntegrationMethod - sample_points::TN - weights::TW -end - -function GaussLegendre(nsamples=20; silence_warnings=true) - if silence_warnings && Bool((nsamples) % 2) - @warn "Must have an even number of GaussLegendre sample points if using for panel self influence" - end - - nodes, weights = FastGaussQuadrature.gausslegendre(nsamples) - - return GaussLegendre(linear_transform((-1, 1), (0, 1), nodes), weights ./ 2.0) -end - -@kwdef struct IntegrationOptions{TN,TS} - nominal::TN = GaussLegendre(20) - singular::TS = GaussLegendre(20) -end diff --git a/src/utilities/types.jl b/src/utilities/types.jl index 7eb70fb0..3263ba17 100644 --- a/src/utilities/types.jl +++ b/src/utilities/types.jl @@ -268,3 +268,43 @@ function quicksolve_options(; ) end +#---------------------------------# +# QUADRATURE TYPES # +#---------------------------------# + +abstract type IntegrationMethod end + +@kwdef struct Romberg{TI,TF,TP,TV} <: IntegrationMethod + max_subdivisions::TI = 10 + atol::TF = 1e-6 +end + +function set_romberg_options(; max_subdivisions=10, atol=1e-6) + return Romberg(; max_subdivisions, atol) +end + +@kwdef struct GaussKronrod{TI,TF} <: IntegrationMethod + order::TI = 5 + maxevals::TI = 1000 + atol::TF = 1e-12 +end + +struct GaussLegendre{TN,TW} <: IntegrationMethod + sample_points::TN + weights::TW +end + +function GaussLegendre(nsamples=20; silence_warnings=true) + if silence_warnings && Bool((nsamples) % 2) + @warn "Must have an even number of GaussLegendre sample points if using for panel self influence" + end + + nodes, weights = FastGaussQuadrature.gausslegendre(nsamples) + + return GaussLegendre(linear_transform((-1, 1), (0, 1), nodes), weights ./ 2.0) +end + +@kwdef struct IntegrationOptions{TN,TS} + nominal::TN = GaussLegendre(20) + singular::TS = GaussLegendre(20) +end From 8e7b040d93bcce7cb7822b344116ef30e62a3b42 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Wed, 3 Apr 2024 19:30:51 -0600 Subject: [PATCH 08/12] remove empty file --- src/types.jl | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 src/types.jl diff --git a/src/types.jl b/src/types.jl deleted file mode 100644 index e69de29b..00000000 From fe9593fe71d69d42839fb6ce545f1d03b1c89287 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Wed, 3 Apr 2024 19:33:15 -0600 Subject: [PATCH 09/12] clean out unused file --- src/preprocess/velocities/integrals.jl | 185 ------------------------- 1 file changed, 185 deletions(-) delete mode 100644 src/preprocess/velocities/integrals.jl diff --git a/src/preprocess/velocities/integrals.jl b/src/preprocess/velocities/integrals.jl deleted file mode 100644 index 64bb010b..00000000 --- a/src/preprocess/velocities/integrals.jl +++ /dev/null @@ -1,185 +0,0 @@ -####################################################################### -## # -## Linear Vortex Panel Integration # -## # -####################################################################### - -##---------------------------------# -## Nominal # -##---------------------------------# - -#""" -#`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -#""" -#function nominal_vortex_panel_integration( -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=[0.0; 1.0], -# debug=false, -#) - -# # Define function to integrate -# function fsample(t) -# return nominal_vortex_induced_velocity_sample( -# t, -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=nondimrange, -# ) -# end - -# V, err = quadgk(fsample, 0.0, 1.0; order=3, maxevals=1e2, atol=1e-6) - -# if debug -# return reshape(V, (2, 2)), err -# # return V, err -# else -# return reshape(V, (2, 2)) -# # return V -# end -#end - -##---------------------------------# -## Singular # -##---------------------------------# - -#""" -#`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -#""" -#function self_vortex_panel_integration( -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=[0.0; 1.0], -# debug=false, -#) - -# # Define function to integrate -# function fsample(t) -# return self_vortex_induced_velocity_sample( -# t, -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=nondimrange, -# ) -# end - -# V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) - -# cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( -# controlpoint[2], influence_length -# ) - -# V .*= influence_length -# V[1:2] .+= cache_vec[1] / 2.0 - -# if debug -# return reshape(V, (2, 2)), err -# else -# return reshape(V, (2, 2)) -# end -#end - -####################################################################### -## # -## Linear Source Panel Integration # -## # -####################################################################### - -##---------------------------------# -## Nominal # -##---------------------------------# - -#""" -#`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -#""" -#function nominal_source_panel_integration( -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=[0.0; 1.0], -# debug=false, -#) - -# # Define function to integrate -# function fsample(t) -# return nominal_source_induced_velocity_sample( -# t, -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=nondimrange, -# ) -# end - -# V, err = quadgk(fsample, 0.0, 1.0; order=3, maxevals=1e2, atol=1e-6) - -# if debug -# return reshape(V, (2, 2)), err -# # return V, err -# else -# return reshape(V, (2, 2)) -# # return V -# end -#end - -##---------------------------------# -## Singular # -##---------------------------------# - -#""" -#`V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] -#""" -#function self_source_panel_integration( -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=[0.0; 1.0], -# debug=false, -#) - -# # Define function to integrate -# function fsample(t) -# return self_source_induced_velocity_sample( -# t, -# node1, -# node2, -# influence_length, -# controlpoint, -# cache_vec; -# nondimrange=nondimrange, -# ) -# end - -# V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) - -# cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( -# controlpoint[2], influence_length -# ) - -# V .*= influence_length -# V[3:4] .+= cache_vec[2] / 2.0 - -# if debug -# return reshape(V, (2, 2)), err -# else -# return reshape(V, (2, 2)) -# end -#end From c97f77225602c474387404ff8f1a8b0e29893893 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Wed, 3 Apr 2024 20:46:25 -0600 Subject: [PATCH 10/12] get integration options in place and ready for unit testing. --- src/DuctAPE.jl | 5 +- src/analysis/analyses.jl | 3 +- src/preprocess/preprocess.jl | 89 ++++++--- src/preprocess/velocities/body_aic.jl | 169 ++++++++++++++--- .../velocities/gausskronrod_integrals.jl | 24 +-- .../velocities/induced_velocity_matrices.jl | 171 ++++++++++++++---- .../velocities/out_of_place_integrals.jl | 102 +++++++++++ .../velocities/romberg_integrals.jl | 31 +--- src/utilities/caches.jl | 29 +++ src/utilities/types.jl | 88 ++++----- 10 files changed, 548 insertions(+), 163 deletions(-) create mode 100644 src/preprocess/velocities/out_of_place_integrals.jl diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index 3b62e3d0..82db5319 100644 --- a/src/DuctAPE.jl +++ b/src/DuctAPE.jl @@ -95,10 +95,11 @@ include("preprocess/velocities/unit_induced_velocities.jl") include("preprocess/velocities/induced_velocity_matrices.jl") include("preprocess/velocities/body_aic.jl") # Quadrature -include("preprocess/velocities/integrals.jl") +include("preprocess/velocities/integrands.jl") +include("preprocess/velocities/out_of_place_integrals.jl") include("preprocess/velocities/gausslegendre_integrals.jl") +include("preprocess/velocities/romberg_integrals.jl") include("preprocess/velocities/gausskronrod_integrals.jl") -include("preprocess/velocities/integrands.jl") ##### ----- PROCESS ----- ##### # Solve and Residual Functions diff --git a/src/analysis/analyses.jl b/src/analysis/analyses.jl index 2923873e..53187c64 100644 --- a/src/analysis/analyses.jl +++ b/src/analysis/analyses.jl @@ -84,7 +84,7 @@ function analyze( # out-of-place version currently has 22,292,181 allocations. # TODO: do this in place for the solve input cache items. eventually will want to have a post-processing and output cache too. - ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps, panels, problem_dimensions = precompute_parameters_iad!( + ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps, panels, problem_dimensions = precompute_parameters!( solve_parameter_tuple.ivr, solve_parameter_tuple.ivw, solve_parameter_tuple.blade_elements, @@ -92,6 +92,7 @@ function analyze( solve_parameter_tuple.wakeK, propulsor; grid_solver_options=options.wake_solver_options, + integration_options=options.integration_options, autoshiftduct=options.autoshiftduct, itcpshift=options.itcpshift, axistol=options.axistol, diff --git a/src/preprocess/preprocess.jl b/src/preprocess/preprocess.jl index 858c5bd3..5e7d86fa 100644 --- a/src/preprocess/preprocess.jl +++ b/src/preprocess/preprocess.jl @@ -242,7 +242,7 @@ end """ """ -function calculate_unit_induced_velocities(problem_dimensions, panels) +function calculate_unit_induced_velocities(problem_dimensions, panels, integration_options) (; nrotor, # number of rotors nwn, # number of wake nodes @@ -277,12 +277,12 @@ function calculate_unit_induced_velocities(problem_dimensions, panels) v_bb=zeros(TF, nbp, nbn, 2), ) - return calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) + return calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_options) end """ """ -function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) +function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_options) # - Extract Tuples - # # Extract induced velocities on rotor (; v_rr, v_rw, v_rb) = ivr @@ -308,6 +308,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) body_vortex_panels.nodemap, body_vortex_panels.influence_length, ones(TF, 2, body_vortex_panels.totpanel), + integration_options, ) # Add influence of body trailing edge gap "panels" @@ -319,6 +320,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) body_vortex_panels.tendotn, body_vortex_panels.tencrossn, body_vortex_panels.teadjnodeidxs, + integration_options, ) # - Rotors on Bodies - # @@ -330,6 +332,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) rotor_source_panels.nodemap, rotor_source_panels.influence_length, ones(TF, 2, rotor_source_panels.totnode), + integration_options, ) # - Wake on Bodies - # @@ -341,6 +344,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) wake_vortex_panels.nodemap, wake_vortex_panels.influence_length, ones(TF, 2, wake_vortex_panels.totpanel), + integration_options, ) # wake "TE panels" to body panels @@ -351,7 +355,8 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) wake_vortex_panels.teinfluence_length, wake_vortex_panels.tendotn, wake_vortex_panels.tencrossn, - wake_vortex_panels.teadjnodeidxs; + wake_vortex_panels.teadjnodeidxs, + integration_options; wake=true, ) @@ -364,6 +369,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) rotor_source_panels.nodemap, rotor_source_panels.influence_length, ones(TF, 2, rotor_source_panels.totnode), + integration_options, ) # - Bodies on Rotors - # @@ -375,6 +381,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) body_vortex_panels.nodemap, body_vortex_panels.influence_length, ones(TF, 2, body_vortex_panels.totpanel), + integration_options, ) # add influence from body trailing edge gap "panels" @@ -383,10 +390,10 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) rotor_source_panels.controlpoint, body_vortex_panels.tenode, body_vortex_panels.teinfluence_length, - # -body_vortex_panels.teinfluence_length, body_vortex_panels.tendotn, body_vortex_panels.tencrossn, body_vortex_panels.teadjnodeidxs, + integration_options, ) # - Wake on Rotors - # @@ -398,6 +405,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) wake_vortex_panels.nodemap, wake_vortex_panels.influence_length, ones(TF, 2, wake_vortex_panels.totpanel), + integration_options, ) # add influence from wake "trailing edge panels" @@ -408,7 +416,8 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) wake_vortex_panels.teinfluence_length, wake_vortex_panels.tendotn, wake_vortex_panels.tencrossn, - wake_vortex_panels.teadjnodeidxs; + wake_vortex_panels.teadjnodeidxs, + integration_options; wake=true, ) @@ -422,6 +431,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) body_vortex_panels.nodemap, body_vortex_panels.influence_length, ones(TF, 2, body_vortex_panels.totpanel), + integration_options, ) # add influence from body trailing edge gap "panels" @@ -434,6 +444,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) body_vortex_panels.tendotn, body_vortex_panels.tencrossn, body_vortex_panels.teadjnodeidxs, + integration_options, ) # - Rotors to Wakes - # @@ -444,6 +455,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) rotor_source_panels.nodemap, rotor_source_panels.influence_length, ones(TF, 2, rotor_source_panels.totpanel), + integration_options, ) # - Wake on Wake - # @@ -455,6 +467,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) wake_vortex_panels.nodemap, wake_vortex_panels.influence_length, ones(TF, 2, wake_vortex_panels.totpanel), + integration_options, ) # add influence from wake "trailing edge panels" @@ -465,7 +478,8 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) wake_vortex_panels.teinfluence_length, wake_vortex_panels.tendotn, wake_vortex_panels.tencrossn, - wake_vortex_panels.teadjnodeidxs; + wake_vortex_panels.teadjnodeidxs, + integration_options; wake=true, ) @@ -475,7 +489,12 @@ end """ """ function initialize_linear_system( - ivb, body_vortex_panels, rotor_source_panels, wake_vortex_panels, Vinf + ivb, + body_vortex_panels, + rotor_source_panels, + wake_vortex_panels, + Vinf, + integration_options, ) # velocities on body @@ -509,7 +528,7 @@ function initialize_linear_system( # Boundary on internal psuedo control point influence coefficients AICpcp = vortex_aic_boundary_on_field( - itcontrolpoint, itnormal, node, nodemap, influence_length + itcontrolpoint, itnormal, node, nodemap, influence_length, integration_options ) # # Add Trailing Edge Gap Panel Influences to panels @@ -534,6 +553,7 @@ function initialize_linear_system( tendotn, tencrossn, teadjnodeidxs, + integration_options, ) # Assemble Raw LHS Matrix into A_bb @@ -565,6 +585,7 @@ function initialize_linear_system( rotor_source_panels.node, rotor_source_panels.nodemap, rotor_source_panels.influence_length, + integration_options, ) # A_pr = calculate_normal_velocity(v_pr, itnormal) @@ -593,6 +614,7 @@ function initialize_linear_system( wake_vortex_panels.node, wake_vortex_panels.nodemap, wake_vortex_panels.influence_length, + integration_options, ) # A_pw = calculate_normal_velocity(v_pw, itnormal) @@ -606,7 +628,8 @@ function initialize_linear_system( wake_vortex_panels.teinfluence_length, wake_vortex_panels.tendotn, wake_vortex_panels.tencrossn, - wake_vortex_panels.teadjnodeidxs; + wake_vortex_panels.teadjnodeidxs, + integration_options; wake=true, ) @@ -623,6 +646,7 @@ function initialize_linear_system!( wake_vortex_panels, Vinf, intermediate_containers, + integration_options, ) # - Clear Containers - # @@ -669,7 +693,13 @@ function initialize_linear_system!( # Boundary on internal psuedo control point influence coefficients vortex_aic_boundary_on_field!( - AICpcp, itcontrolpoint, itnormal, node, nodemap, influence_length + AICpcp, + itcontrolpoint, + itnormal, + node, + nodemap, + influence_length, + integration_options, ) # Add Trailing Edge Gap Panel Influences to internal pseudo control point @@ -682,6 +712,7 @@ function initialize_linear_system!( tendotn, tencrossn, teadjnodeidxs, + integration_options, ) # Assemble Raw LHS Matrix into A_bb @@ -720,6 +751,7 @@ function initialize_linear_system!( rotor_source_panels.node, rotor_source_panels.nodemap, rotor_source_panels.influence_length, + integration_options, ) ##### ----- Wake AIC ----- ##### @@ -734,6 +766,7 @@ function initialize_linear_system!( wake_vortex_panels.node, wake_vortex_panels.nodemap, wake_vortex_panels.influence_length, + integration_options, ) # add contributions from wake "trailing edge panels" on pseudo control point @@ -745,7 +778,8 @@ function initialize_linear_system!( wake_vortex_panels.teinfluence_length, wake_vortex_panels.tendotn, wake_vortex_panels.tencrossn, - wake_vortex_panels.teadjnodeidxs; + wake_vortex_panels.teadjnodeidxs, + integration_options; wake=true, ) @@ -1267,6 +1301,7 @@ end function precompute_parameters( propulsor; grid_solver_options=GridSolverOptions(), + integration_options=IntegrationOptions(), autoshiftduct=true, itcpshift=0.05, axistol=1e-15, @@ -1321,7 +1356,6 @@ function precompute_parameters( ) return precompute_parameters( - problem_dimensions, rp_duct_coordinates, rp_centerbody_coordinates, wake_grid, @@ -1331,7 +1365,9 @@ function precompute_parameters( rotorstator_parameters, paneling_constants, operating_point, - reference_parameters; + reference_parameters, + integration_options, + problem_dimensions; itcpshift=itcpshift, axistol=axistol, tegaptol=tegaptol, @@ -1354,6 +1390,7 @@ function precompute_parameters( paneling_constants, operating_point, reference_parameters, + integration_options, problem_dimensions=nothing; autoshiftduct=true, itcpshift=0.05, @@ -1386,7 +1423,9 @@ function precompute_parameters( # - Compute Influence Matrices - # ivr, ivw, ivb = calculate_unit_induced_velocities( - problem_dimensions, (; body_vortex_panels, rotor_source_panels, wake_vortex_panels) + problem_dimensions, + (; body_vortex_panels, rotor_source_panels, wake_vortex_panels), + integration_options, ) # - Set up Linear System - # @@ -1396,6 +1435,7 @@ function precompute_parameters( rotor_source_panels, wake_vortex_panels, operating_point.Vinf, + integration_options, ) # - Interpolate Blade Elements - # @@ -1448,8 +1488,8 @@ function precompute_parameters!( linsys, wakeK, propulsor; - #TODO: put in the actual defaults here grid_solver_options=GridSolverOptions(), + integration_options=IntegrationOptions(), autoshiftduct=true, itcpshift=0.05, axistol=1e-15, @@ -1512,6 +1552,7 @@ function precompute_parameters!( reference_parameters, problem_dimensions; grid_solver_options=grid_solver_options, + integration_options=integration_options, autoshiftduct=autoshiftduct, itcpshift=itcpshift, axistol=axistol, @@ -1540,6 +1581,7 @@ function precompute_parameters!( reference_parameters, problem_dimensions=nothing; grid_solver_options=options.wake_solver_options, + integration_options=options.integration_options, autoshiftduct=options.autoshiftduct, itcpshift=options.itcpshift, axistol=options.axistol, @@ -1607,7 +1649,11 @@ function precompute_parameters!( v_bw=zeros(TF, problem_dimensions.nbp, problem_dimensions.nwn, 2), ) calculate_unit_induced_velocities!( - ivr, ivw, ivb, (; body_vortex_panels, rotor_source_panels, wake_vortex_panels) + ivr, + ivw, + ivb, + (; body_vortex_panels, rotor_source_panels, wake_vortex_panels), + integration_options, ) # - Set up Linear System - # @@ -1628,6 +1674,7 @@ function precompute_parameters!( wake_vortex_panels, operating_point.Vinf[1], intermediate_containers, + integration_options ) # - Interpolate Blade Elements - # @@ -1734,12 +1781,12 @@ end # # - Compute Influence Matrices - # # # TODO: test this function -# calculate_unit_induced_velocities!(ivr, ivw, ivb, panels) +# calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_options) # # - Set up Linear System - # # # TODO: test this function # A_bb_LU, lu_decomp_flag = initialize_linear_system!( -# linsys, ivb, panels.body_vortex_panels, AICn, AICpcp +# linsys, ivb, panels.body_vortex_panels, AICn, AICpcp, integration_options # ) # # - Interpolate Blade Elements - # @@ -2049,7 +2096,7 @@ function initialize_strengths!( # gamb = ImplicitAD.implicit_linear( # linsys.A_bb, copy(linsys.b_bf); lsolve=ldiv!, Af=linsys.A_bb_LU # ) - gamb = zeros(size(ivr.v_rb,2)+2) + gamb = zeros(size(ivr.v_rb, 2) + 2) # - Get body-induced velocities on rotors - # vzb = zeros(TF, nbe, nrotor) @@ -2207,7 +2254,7 @@ function initialize_strengths!( # Gamr struggles to converge if it's not initially positive... for g in eachindex(Gamr) - if Gamr[g]<0 + if Gamr[g] < 0 Gamr[g] = 0.05 end end diff --git a/src/preprocess/velocities/body_aic.jl b/src/preprocess/velocities/body_aic.jl index 5bdfa8d8..73091a11 100644 --- a/src/preprocess/velocities/body_aic.jl +++ b/src/preprocess/velocities/body_aic.jl @@ -24,7 +24,7 @@ Used for constructing the LHS influence Matrix for the panel method system, as w - `AICn::Matrix{Float}` : N controlpoint x N+1 node array of V dot nhat values """ function vortex_aic_boundary_on_boundary( - controlpoint, normal, node, nodemap, influence_length + controlpoint, normal, node, nodemap, influence_length, integration_options ) T = promote_type(eltype(node), eltype(controlpoint)) M = size(controlpoint, 2) @@ -33,7 +33,7 @@ function vortex_aic_boundary_on_boundary( AICn = zeros(T, M, N) vortex_aic_boundary_on_boundary!( - AICn, controlpoint, normal, node, nodemap, influence_length + AICn, controlpoint, normal, node, nodemap, influence_length, integration_options ) return AICn @@ -57,13 +57,29 @@ Used for constructing the LHS influence Matrix for the panel method system, as w - `influence_length::Vector{Float}` : lengths over which vortex ring influence is applied on the surface. """ function vortex_aic_boundary_on_boundary!( - AICn, controlpoint, normal, node, nodemap, influence_length; cache_vec=nothing + AICn, + controlpoint, + normal, + node, + nodemap, + influence_length, + integration_options; + integration_caches=nothing, ) # NOTE: it is slighlty faster/fewer allocations to just define a new static array in the loop than to preallocate such a small matrix. # vel = zeros(eltype(AICn), 2, 2) - if isnothing(cache_vec) - cache_vec = zeros(eltype(node), 20) + if isnothing(integration_caches) + # integration_cache = zeros(eltype(controlpoint), 20) + nominal_integration_cache = allocate_integration_containers( + integration_options.nominal, eltype(VEL) + ) + singular_integration_cache = allocate_integration_containers( + integration_options.singular, eltype(VEL) + ) + else + nominal_integration_cache = integration_caches.nominal + singular_integration_cache = integration_caches.singular end # loop through panels doing the influencing @@ -77,12 +93,26 @@ function vortex_aic_boundary_on_boundary!( if i != j # vel .= nominal_vortex_panel_integration(n1, n2, lj, cpi) vel = StaticArrays.SMatrix{2,2}( - nominal_vortex_panel_integration(n1, n2, lj, cpi, cache_vec) + nominal_vortex_panel_integration( + integration_options.nominal, + n1, + n2, + lj, + cpi, + nominal_integration_cache, + ), ) else # vel .= self_vortex_panel_integration(n1, n2, lj, cpi) vel = StaticArrays.SMatrix{2,2}( - self_vortex_panel_integration(n1, n2, lj, cpi, cache_vec) + self_vortex_panel_integration( + integration_options.singular, + n1, + n2, + lj, + cpi, + singular_integration_cache, + ), ) end @@ -112,7 +142,13 @@ Used for constructing the LHS influence Matrix for the panel method system, as w - `AICn::Matrix{Float}` : N controlpoint x N+1 node array of V dot nhat values """ function vortex_aic_boundary_on_field( - controlpoint, normal, node, nodemap, influence_length; cache_vec=nothing + controlpoint, + normal, + node, + nodemap, + influence_length, + integration_options; + integration_caches=nothing, ) T = promote_type(eltype(node), eltype(controlpoint)) M = size(controlpoint, 2) @@ -121,7 +157,14 @@ function vortex_aic_boundary_on_field( AICn = zeros(T, M, N) vortex_aic_boundary_on_field!( - AICn, controlpoint, normal, node, nodemap, influence_length; cache_vec=cache_vec + AICn, + controlpoint, + normal, + node, + nodemap, + influence_length, + integration_options; + integration_caches=integration_cache, ) return AICn @@ -140,13 +183,28 @@ Used for constructing the LHS influence Matrix for the panel method system, as w - `influence_length::Vector{Float}` : lengths over which vortex ring influence is applied on the surface. """ function vortex_aic_boundary_on_field!( - AICn, controlpoint, normal, node, nodemap, influence_length; cache_vec=nothing + AICn, + controlpoint, + normal, + node, + nodemap, + influence_length, + integration_options; + integration_caches=nothing, ) # vel = zeros(eltype(AICn), 2, 2) - if isnothing(cache_vec) - cache_vec = zeros(eltype(node), 20) + if isnothing(integration_caches) + # integration_cache = zeros(eltype(controlpoint), 20) + nominal_integration_cache = allocate_integration_containers( + integration_options.nominal, eltype(AICn) + ) + singular_integration_cache = allocate_integration_containers( + integration_options.singular, eltype(AICn) + ) + else + nominal_integration_cache = integration_caches.nominal + singular_integration_cache = integration_caches.singular end - # Loop through control points being influenced for (i, (cpi, nhat)) in enumerate(zip(eachcol(controlpoint), eachcol(normal))) # loop through panels doing the influencing @@ -159,13 +217,27 @@ function vortex_aic_boundary_on_field!( # if so: # vel .= self_vortex_panel_integration(n1, n2, lj, cpi) vel = StaticArrays.SMatrix{2,2}( - self_vortex_panel_integration(n1, n2, lj, cpi, cache_vec) + self_vortex_panel_integration( + integration_options.singular, + n1, + n2, + lj, + cpi, + singular_integration_cache, + ), ) else # if not: # vel .= nominal_vortex_panel_integration(n1, n2, lj, cpi) vel = StaticArrays.SMatrix{2,2}( - nominal_vortex_panel_integration(n1, n2, lj, cpi, cache_vec) + nominal_vortex_panel_integration( + integration_options.nominal, + n1, + n2, + lj, + cpi, + nominal_integration_cache, + ), ) end @@ -221,12 +293,18 @@ function add_te_gap_aic!( teinfluence_length, tendotn, tencrossn, - teadjnodeidxs; + teadjnodeidxs, + integration_options; wake=false, - cache_vec=nothing, + integration_caches=nothing, ) - if isnothing(cache_vec) - cache_vec = zeros(eltype(controlpoint), 20) + if isnothing(integration_caches) + # integration_cache = zeros(eltype(controlpoint), 20) + nominal_integration_cache = allocate_integration_containers( + integration_options.nominal, eltype(AICn) + ) + else + nominal_integration_cache = integration_caches end # Loop through control points being influenced @@ -244,12 +322,22 @@ function add_te_gap_aic!( # get unit induced velocity from the panel onto the control point vvel = StaticArrays.SMatrix{2,2}( nominal_vortex_panel_integration( - tenode[j, 1, :], tenode[j, 2, :], lj, cpi, cache_vec + integration_options.nominal, + tenode[j, 1, :], + tenode[j, 2, :], + lj, + cpi, + nominal_integration_cache, ), ) svel = StaticArrays.SMatrix{2,2}( nominal_source_panel_integration( - tenode[j, 1, :], tenode[j, 2, :], lj, cpi, cache_vec + integration_options.nominal, + tenode[j, 1, :], + tenode[j, 2, :], + lj, + cpi, + nominal_integration_cache, ), ) @@ -286,7 +374,13 @@ Used for constructing the RHS boundary conditions due to rotor source panels. - `AICn::Matrix{Float}` : N controlpoint x N+1 node array of V dot nhat values """ function source_aic( - controlpoint, normal, node, nodemap, influence_length; cache_vec=nothing + controlpoint, + normal, + node, + nodemap, + influence_length, + integration_options; + integration_caches=nothing, ) T = promote_type(eltype(node), eltype(controlpoint)) M = size(controlpoint, 2) @@ -295,7 +389,14 @@ function source_aic( AICn = zeros(T, M, N) source_aic!( - AICn, controlpoint, normal, node, nodemap, influence_length; cache_vec=cache_vec + AICn, + controlpoint, + normal, + node, + nodemap, + influence_length, + integration_options; + integration_caches=integration_caches, ) return AICn @@ -315,11 +416,23 @@ end - `influence_length::Vector{Float}` : lengths over which vortex ring influence is applied on the surface. """ function source_aic!( - AICn, controlpoint, normal, node, nodemap, influence_length; cache_vec=nothing + AICn, + controlpoint, + normal, + node, + nodemap, + influence_length, + integration_options; + integration_caches=nothing, ) # vel = zeros(eltype(AICn), 2, 2) - if isnothing(cache_vec) - cache_vec = zeros(eltype(node), 20) + if isnothing(integration_caches) + # integration_cache = zeros(eltype(controlpoint), 20) + nominal_integration_cache = allocate_integration_containers( + integration_options.nominal, eltype(AICn) + ) + else + nominal_integration_cache = integration_caches end # Loop through control points being influenced @@ -332,7 +445,9 @@ function source_aic!( # get unit induced velocity from the panel onto the control point # vel .= nominal_vortex_panel_integration(n1, n2, lj, cpi) vel = StaticArrays.SMatrix{2,2}( - nominal_source_panel_integration(n1, n2, lj, cpi, cache_vec) + nominal_source_panel_integration( + integration_options.nominal, n1, n2, lj, cpi, nominal_integration_cache + ), ) for k in 1:2 diff --git a/src/preprocess/velocities/gausskronrod_integrals.jl b/src/preprocess/velocities/gausskronrod_integrals.jl index e76f3090..d547ff70 100644 --- a/src/preprocess/velocities/gausskronrod_integrals.jl +++ b/src/preprocess/velocities/gausskronrod_integrals.jl @@ -11,14 +11,14 @@ function nominal_vortex_panel_integration( node2, influence_length, controlpoint, - cache_vec; + sample_cache; debug=false, ) # Define function to integrate function fsample(t) return nominal_vortex_induced_velocity_sample( - t, node1, node2, influence_length, controlpoint, cache_vec + t, node1, node2, influence_length, controlpoint, sample_cache ) end @@ -49,14 +49,14 @@ function self_vortex_panel_integration( node2, influence_length, controlpoint, - cache_vec; + sample_cache; debug=false, ) # Define function to integrate function fsample(t) return self_vortex_induced_velocity_sample( - t, node1, node2, influence_length, controlpoint, cache_vec; + t, node1, node2, influence_length, controlpoint, sample_cache; ) end @@ -70,12 +70,12 @@ function self_vortex_panel_integration( atol=integration_options.atol, ) - cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( + sample_cache[1], sample_cache[2] = analytically_integrated_vortex_influence( controlpoint[2], influence_length ) V .*= influence_length - V[1:2] .+= cache_vec[1] / 2.0 + V[1:2] .+= sample_cache[1] / 2.0 if debug return reshape(V, (2, 2)), err @@ -97,14 +97,14 @@ function nominal_source_panel_integration( node2, influence_length, controlpoint, - cache_vec; + sample_cache; debug=false, ) # Define function to integrate function fsample(t) return nominal_source_induced_velocity_sample( - t, node1, node2, influence_length, controlpoint, cache_vec; + t, node1, node2, influence_length, controlpoint, sample_cache; ) end @@ -135,14 +135,14 @@ function self_source_panel_integration( node2, influence_length, controlpoint, - cache_vec; + sample_cache; debug=false, ) # Define function to integrate function fsample(t) return self_source_induced_velocity_sample( - t, node1, node2, influence_length, controlpoint, cache_vec; + t, node1, node2, influence_length, controlpoint, sample_cache; ) end @@ -156,12 +156,12 @@ function self_source_panel_integration( atol=integration_options.atol, ) - cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( + sample_cache[1], sample_cache[2] = analytically_integrated_source_influence( controlpoint[2], influence_length ) V .*= influence_length - V[3:4] .+= cache_vec[2] / 2.0 + V[3:4] .+= sample_cache[2] / 2.0 if debug return reshape(V, (2, 2)), err diff --git a/src/preprocess/velocities/induced_velocity_matrices.jl b/src/preprocess/velocities/induced_velocity_matrices.jl index 949c3d17..d575b6fc 100644 --- a/src/preprocess/velocities/induced_velocity_matrices.jl +++ b/src/preprocess/velocities/induced_velocity_matrices.jl @@ -1,8 +1,3 @@ -#= -Influence matrices for panels on points -given in terms of raw velocity components -=# - ###################################################################### # # # Matrix of velocities induced by panels on points # @@ -28,7 +23,13 @@ Used for getting the unit induced velocities due to the body panels on the rotor - `AIC::Array{Float}` : N-controlpoint x N-node x [vz, vr] array of induced velocity components """ function induced_velocities_from_vortex_panels_on_points( - controlpoints, nodes, nodemap, influence_lengths, strengths; cache_vec=nothing + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + integration_options; + integration_caches=nothing, ) # Initialize @@ -41,7 +42,8 @@ function induced_velocities_from_vortex_panels_on_points( nodes, nodemap, influence_lengths, - strengths; + strengths, + integration_options; cache_vec=cache_vec, ) @@ -61,11 +63,27 @@ Used for getting the unit induced velocities due to the body panels on the rotor - `gamma::Vector{Float}` : vortex constant circulation values """ function induced_velocities_from_vortex_panels_on_points!( - VEL, controlpoint, node, nodemap, influence_length, strength; cache_vec=nothing + VEL, + controlpoint, + node, + nodemap, + influence_length, + strength, + integration_options; + integration_caches=nothing, ) - # vel = zeros(eltype(VEL), 2, 2) - if isnothing(cache_vec) - cache_vec = zeros(eltype(node), 20) + # vel = zeros(VEL, 2, 2) + if isnothing(integration_caches) + #integration_cache = zeros(eltype(node), 20) + nominal_integration_cache = allocate_integration_containers( + integration_options.nominal, VEL + ) + singular_integration_cache = allocate_integration_containers( + integration_options.singular, VEL + ) + else + nominal_integration_cache = integration_caches.nominal + singular_integration_cache = integration_caches.singular end # loop through panels doing the influencing @@ -81,13 +99,27 @@ function induced_velocities_from_vortex_panels_on_points!( # if so: # vel .= self_vortex_panel_integration(n1, n2, lj, cpi) vel = StaticArrays.SMatrix{2,2}( - self_vortex_panel_integration(n1, n2, lj, cpi, cache_vec) + self_vortex_panel_integration( + integration_options.singular, + n1, + n2, + lj, + cpi, + singular_integration_cache, + ), ) else # if not: # vel .= nominal_vortex_panel_integration(n1, n2, lj, cpi) vel = StaticArrays.SMatrix{2,2}( - nominal_vortex_panel_integration(n1, n2, lj, cpi, cache_vec) + nominal_vortex_panel_integration( + integration_options.nominal, + n1, + n2, + lj, + cpi, + nominal_integration_cache, + ), ) end @@ -120,7 +152,13 @@ Used for getting the unit induced velocities due to the body panels on the rotor - `AIC::Array{Float}` : N-controlpoint x N-node x [vz, vr] array of induced velocity components """ function induced_velocities_from_source_panels_on_points( - controlpoints, nodes, nodemap, influence_lengths, strengths; cache_vec=nothing + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + integration_options; + integration_caches=nothing, ) # Initialize @@ -132,9 +170,10 @@ function induced_velocities_from_source_panels_on_points( controlpoints, nodes, nodemap, - influence_lengths, - strengths; - cache_vec=cache_vec, + influence_lengthss, + strengths, + integration_options; + integration_caches=integration_caches, ) return VEL @@ -153,15 +192,29 @@ Used for getting the unit induced velocities due to the body panels on the rotor - `gamma::Vector{Float}` : source constant circulation values """ function induced_velocities_from_source_panels_on_points!( - VEL, controlpoint, node, nodemap, influence_length, strength; cache_vec=nothing + VEL, + controlpoint, + node, + nodemap, + influence_length, + strength, + integration_options; + integration_caches=nothing, ) - # vel = zeros(eltype(VEL), 2, 2) - if isnothing(cache_vec) - cache_vec = zeros(eltype(node), 20) + # vel = zeros(VEL, 2, 2) + if isnothing(integration_caches) + #integration_cache = zeros(eltype(node), 20) + nominal_integration_cache = allocate_integration_containers( + integration_options.nominal, VEL + ) + singular_integration_cache = allocate_integration_containers( + integration_options.singular, VEL + ) + else + nominal_integration_cache = integration_caches.nominal + singular_integration_cache = integration_caches.singular end - #TODO; for speedups, update panel initialization to flip rows and columns such that these functions use eachcol rather than eachrow - # loop through panels doing the influencing for (j, (nmap, lj, sigmaj)) in enumerate(zip(eachcol(nodemap), influence_length, eachcol(strength))) @@ -176,13 +229,27 @@ function induced_velocities_from_source_panels_on_points!( # if so: # vel .= self_source_panel_integration(n1, n2, lj, cpi) vel = StaticArrays.SMatrix{2,2}( - self_source_panel_integration(n1, n2, lj, cpi, cache_vec) + self_source_panel_integration( + integration_options.singular, + n1, + n2, + lj, + cpi, + singular_integration_cache, + ), ) else # if not: # vel .= nominal_source_panel_integration(n1, n2, lj, cpi) vel = StaticArrays.SMatrix{2,2}( - nominal_source_panel_integration(n1, n2, lj, cpi, cache_vec) + nominal_source_panel_integration( + integration_options.nominal, + n1, + n2, + lj, + cpi, + nominal_integration_cache, + ), ) end @@ -208,15 +275,25 @@ function induced_velocities_from_trailing_edge_gap_panel!( teinfluence_length, tendotn, tencrossn, - teadjnodeidxs; + teadjnodeidxs, + integration_options; wake=false, - cache_vec=nothing, + integration_caches=nothing, ) # vvel = zeros(eltype(AICn), 2, 2) # svel = zeros(eltype(AICn), 2, 2) - if isnothing(cache_vec) - cache_vec = zeros(eltype(controlpoint), 20) + if isnothing(integration_caches) + # integration_cache = zeros(eltype(controlpoint), 20) + nominal_integration_cache = allocate_integration_containers( + integration_options.nominal, VEL + ) + singular_integration_cache = allocate_integration_containers( + integration_options.singular, VEL + ) + else + nominal_integration_cache = integration_caches.nominal + singular_integration_cache = integration_caches.singular end # Loop through control points being influenced @@ -237,18 +314,46 @@ function induced_velocities_from_trailing_edge_gap_panel!( if isapprox(cpi, 0.5 * (n1 .+ n2)) # if so: vvel = StaticArrays.SMatrix{2,2}( - self_vortex_panel_integration(n1, n2, lj, cpi, cache_vec) + self_vortex_panel_integration( + integration_options.singular, + n1, + n2, + lj, + cpi, + singular_integration_cache, + ), ) svel = StaticArrays.SMatrix{2,2}( - self_source_panel_integration(n1, n2, lj, cpi, cache_vec) + self_source_panel_integration( + integration_options.singular, + n1, + n2, + lj, + cpi, + nominal_integration_cache, + ), ) else # if not: vvel = StaticArrays.SMatrix{2,2}( - nominal_vortex_panel_integration(n1, n2, lj, cpi, cache_vec) + nominal_vortex_panel_integration( + integration_options.nominal, + n1, + n2, + lj, + cpi, + singular_integration_cache, + ), ) svel = StaticArrays.SMatrix{2,2}( - nominal_source_panel_integration(n1, n2, lj, cpi, cache_vec) + nominal_source_panel_integration( + integration_options.nominal, + n1, + n2, + lj, + cpi, + nominal_integration_cache, + ), ) end diff --git a/src/preprocess/velocities/out_of_place_integrals.jl b/src/preprocess/velocities/out_of_place_integrals.jl new file mode 100644 index 00000000..f4b52a0f --- /dev/null +++ b/src/preprocess/velocities/out_of_place_integrals.jl @@ -0,0 +1,102 @@ +#---------------------------------# +# VORTEX # +#---------------------------------# + +""" +""" +function nominal_vortex_panel_integration( + integration_options, + node1, + node2, + influence_length, + controlpoint, + containers; + debug=false, +) + TF = promote_type(eltype(node1), eltype(node2)) + V = zeros(TF, 2, 2) + return nominal_vortex_panel_integration!( + integration_options, + V, + node1, + node2, + influence_length, + controlpoint, + containers; + debug=debug, + ) +end + +function self_vortex_panel_integration( + integration_options, + node1, + node2, + influence_length, + controlpoint, + sample_cache; + debug=false, +) + TF = promote_type(eltype(node1), eltype(node2)) + V = zeros(TF, 2, 2) + return self_vortex_panel_integration!( + integration_options, + V, + node1, + node2, + influence_length, + controlpoint, + sample_cache; + debug=debug, + ) +end + +#---------------------------------# +# SOURCE # +#---------------------------------# + +function nominal_source_panel_integration( + integration_options, + node1, + node2, + influence_length, + controlpoint, + containers; + debug=false, +) + TF = promote_type(eltype(node1), eltype(node2)) + V = zeros(TF, 2, 2) + return nominal_source_panel_integration!( + integration_options, + V, + node1, + node2, + influence_length, + controlpoint, + containers; + debug=debug, + ) +end + +function self_source_panel_integration( + integration_options, + node1, + node2, + influence_length, + controlpoint, + sample_cache; + debug=false, +) + TF = promote_type(eltype(node1), eltype(node2)) + V = zeros(TF, 2, 2) + + return self_source_panel_integration!( + integration_options, + V, + node1, + node2, + influence_length, + controlpoint, + sample_cache; + debug=debug, + ) +end diff --git a/src/preprocess/velocities/romberg_integrals.jl b/src/preprocess/velocities/romberg_integrals.jl index f8dfd4da..db5e63ab 100644 --- a/src/preprocess/velocities/romberg_integrals.jl +++ b/src/preprocess/velocities/romberg_integrals.jl @@ -10,7 +10,7 @@ the type of `h` and `T` is the type of the extrapolated `f(0)` **result**. This should be a floating-point type, i.e. `fh` should contain `float(f(h))` if the function you are extrapolating is not already floating-point-valued.) """ -function extrapolate!(V, err, fh; power=2, atol=1e-6) where {TF1,TF2} +function extrapolate!(V, err, fh; power=2, atol=1e-6) # - rename for convenience - # (f0, h) = first(fh) @@ -53,29 +53,6 @@ function extrapolate!(V, err, fh; power=2, atol=1e-6) where {TF1,TF2} return V, err end -function allocate_integration_containers( - integration_options::Romberg, dispatch_type::TF; cache_size=20 -) where {TF} - return (; - sample_cache=zeros(TF, cache_size), - samples=[(zeros(TF, 4), TF[0.0]) for i in 1:(integration_options.max_subdivisions)], - ) -end - -function allocate_integration_containers( - integration_options::GaussLegendre, dispatch_type::TF; cache_size=20 -) where {TF} - return (; - sample_cache=zeros(TF, cache_size), samples=zeros(TF, integration_options.nsamples) - ) -end - -function allocate_integration_containers( - integration_options::GaussKronrod, dispatch_type::TF; cache_size=20 -) where {TF} - return (; sample_cache=zeros(TF, cache_size)) -end - #---------------------------------# # VORTEX # #---------------------------------# @@ -139,8 +116,9 @@ end """ `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ -function self_vortex_panel_integration( +function self_vortex_panel_integration!( integration_options::Romberg, + V, node1, node2, influence_length, @@ -269,8 +247,9 @@ end """ `V::Matrix{Float}` : velocity components due to the jth and j+1th nodes in the format: [vz_j vr_j; vz_{j+1} vr_{j+1}] """ -function self_source_panel_integration( +function self_source_panel_integration!( integration_options::Romberg, + V, node1, node2, influence_length, diff --git a/src/utilities/caches.jl b/src/utilities/caches.jl index f7339892..252ed03d 100644 --- a/src/utilities/caches.jl +++ b/src/utilities/caches.jl @@ -1358,3 +1358,32 @@ function withdraw_post_parameter_cache(vec, dims) return (;) end +#---------------------------------# +# Integration Containers # +#---------------------------------# +function allocate_integration_containers( + integration_options::Romberg, dispatch_type; cache_size=20 +) + TF = eltype(dispatch_type) + return (; + sample_cache=zeros(TF, cache_size), + samples=[(zeros(TF, 4), TF[0.0]) for i in 1:(integration_options.max_subdivisions)], + ) +end + +function allocate_integration_containers( + integration_options::GaussLegendre, dispatch_type; cache_size=20 +) + TF = eltype(dispatch_type) + return (; + sample_cache=zeros(TF, cache_size), + samples=zeros(TF, length(integration_options.sample_points), 4), + ) +end + +function allocate_integration_containers( + integration_options::GaussKronrod, dispatch_type; cache_size=20 +) + TF = eltype(dispatch_type) + return zeros(TF, cache_size) +end diff --git a/src/utilities/types.jl b/src/utilities/types.jl index 3263ba17..b3cc4038 100644 --- a/src/utilities/types.jl +++ b/src/utilities/types.jl @@ -16,6 +16,9 @@ abstract type ExternalSolverOptions <: SolverOptionsType end # - Wake Solver Options - # abstract type GridSolverOptionsType end +# - Quadrature Types - # +abstract type IntegrationMethod end + #---------------------------------# # INPUT TYPES # #---------------------------------# @@ -149,6 +152,45 @@ function verify_input(propulsor) # TODO: go find all the various asserts and put them here end +#---------------------------------# +# QUADRATURE TYPES # +#---------------------------------# + +@kwdef struct Romberg{TI,TF,TP,TV} <: IntegrationMethod + max_subdivisions::TI = 10 + atol::TF = 1e-6 +end + +function set_romberg_options(; max_subdivisions=10, atol=1e-6) + return Romberg(; max_subdivisions, atol) +end + +@kwdef struct GaussKronrod{TI,TF} <: IntegrationMethod + order::TI = 5 + maxevals::TI = 1000 + atol::TF = 1e-12 +end + +struct GaussLegendre{TN,TW} <: IntegrationMethod + sample_points::TN + weights::TW +end + +function GaussLegendre(nsamples=20; silence_warnings=true) + if silence_warnings && Bool((nsamples) % 2) + @warn "Must have an even number of GaussLegendre sample points if using for panel self influence" + end + + nodes, weights = FastGaussQuadrature.gausslegendre(nsamples) + + return GaussLegendre(linear_transform((-1, 1), (0, 1), nodes), weights ./ 2.0) +end + +@kwdef struct IntegrationOptions{TN,TS} + nominal::TN = GaussLegendre(20) + singular::TS = GaussLegendre(20) +end + #---------------------------------# # OPTION TYPES # #---------------------------------# @@ -226,7 +268,9 @@ end # - Full Option Set - # """ """ -@kwdef struct Options{TB,TF,TS,Tin,TSo<:SolverOptionsType,WS<:GridSolverOptionsType} +@kwdef struct Options{ + TB,TF,TS,Tin,TIo<:IntegrationOptions,TSo<:SolverOptionsType,WS<:GridSolverOptionsType +} # - General Options - # verbose::TB = false silence_warnings::TB = true @@ -238,6 +282,8 @@ end itcpshift::TF = 0.05 axistol::TF = 1e-15 tegaptol::TF = 1e1 * eps() + # - Integration Options - # + integration_options::TIo = IntegrationOptions() # - Post-processing Options - # write_outputs::TB = false outfile::TS = "outputs.jl" @@ -268,43 +314,3 @@ function quicksolve_options(; ) end -#---------------------------------# -# QUADRATURE TYPES # -#---------------------------------# - -abstract type IntegrationMethod end - -@kwdef struct Romberg{TI,TF,TP,TV} <: IntegrationMethod - max_subdivisions::TI = 10 - atol::TF = 1e-6 -end - -function set_romberg_options(; max_subdivisions=10, atol=1e-6) - return Romberg(; max_subdivisions, atol) -end - -@kwdef struct GaussKronrod{TI,TF} <: IntegrationMethod - order::TI = 5 - maxevals::TI = 1000 - atol::TF = 1e-12 -end - -struct GaussLegendre{TN,TW} <: IntegrationMethod - sample_points::TN - weights::TW -end - -function GaussLegendre(nsamples=20; silence_warnings=true) - if silence_warnings && Bool((nsamples) % 2) - @warn "Must have an even number of GaussLegendre sample points if using for panel self influence" - end - - nodes, weights = FastGaussQuadrature.gausslegendre(nsamples) - - return GaussLegendre(linear_transform((-1, 1), (0, 1), nodes), weights ./ 2.0) -end - -@kwdef struct IntegrationOptions{TN,TS} - nominal::TN = GaussLegendre(20) - singular::TS = GaussLegendre(20) -end From 9a81bfa735345c2627830f53140d13d67a5a4f34 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Thu, 4 Apr 2024 10:30:42 -0600 Subject: [PATCH 11/12] update code and tests to allow for various implemented quadrature options --- src/DuctAPE.jl | 2 - src/preprocess/preprocess.jl | 18 +- src/preprocess/velocities/body_aic.jl | 6 +- .../velocities/gausslegendre_integrals.jl | 5 + .../velocities/induced_velocity_matrices.jl | 4 +- .../velocities/romberg_integrals.jl | 29 +- src/utilities/types.jl | 12 +- src/utilities/utils.jl | 12 +- test/induced_velocities.jl | 770 +++++++++++++++++- test/linear_system_assembly.jl | 7 + test/runtests.jl | 5 +- 11 files changed, 795 insertions(+), 75 deletions(-) diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index 82db5319..989e5ec8 100644 --- a/src/DuctAPE.jl +++ b/src/DuctAPE.jl @@ -41,8 +41,6 @@ using ForwardDiff # used for jacobian for newton solver using FLOWMath # used for various items, mostly interpolation using Printf # used when verbose option is selected -using Primes # used in Romberg integration settings - #---------------------------------# # EXPORTS # #---------------------------------# diff --git a/src/preprocess/preprocess.jl b/src/preprocess/preprocess.jl index 5e7d86fa..924c692e 100644 --- a/src/preprocess/preprocess.jl +++ b/src/preprocess/preprocess.jl @@ -1580,15 +1580,15 @@ function precompute_parameters!( operating_point, reference_parameters, problem_dimensions=nothing; - grid_solver_options=options.wake_solver_options, - integration_options=options.integration_options, - autoshiftduct=options.autoshiftduct, - itcpshift=options.itcpshift, - axistol=options.axistol, - tegaptol=options.tegaptol, - finterp=options.finterp, - silence_warnings=options.silence_warnings, - verbose=options.verbose, + grid_solver_options=GridSolverOptions(), + integration_options=IntegrationOptions(), + autoshiftduct=true, + itcpshift=0.05, + axistol=1e-15, + tegaptol=1e1 * eps(), + finterp=fm.akima, + silence_warnings=true, + verbose=false, ) # - Reset Caches - # diff --git a/src/preprocess/velocities/body_aic.jl b/src/preprocess/velocities/body_aic.jl index 73091a11..53f90079 100644 --- a/src/preprocess/velocities/body_aic.jl +++ b/src/preprocess/velocities/body_aic.jl @@ -72,10 +72,10 @@ function vortex_aic_boundary_on_boundary!( if isnothing(integration_caches) # integration_cache = zeros(eltype(controlpoint), 20) nominal_integration_cache = allocate_integration_containers( - integration_options.nominal, eltype(VEL) + integration_options.nominal, eltype(AICn) ) singular_integration_cache = allocate_integration_containers( - integration_options.singular, eltype(VEL) + integration_options.singular, eltype(AICn) ) else nominal_integration_cache = integration_caches.nominal @@ -164,7 +164,7 @@ function vortex_aic_boundary_on_field( nodemap, influence_length, integration_options; - integration_caches=integration_cache, + integration_caches=integration_caches, ) return AICn diff --git a/src/preprocess/velocities/gausslegendre_integrals.jl b/src/preprocess/velocities/gausslegendre_integrals.jl index 546e4141..1c1ff6a0 100644 --- a/src/preprocess/velocities/gausslegendre_integrals.jl +++ b/src/preprocess/velocities/gausslegendre_integrals.jl @@ -15,6 +15,8 @@ function nominal_vortex_panel_integration!( containers; debug=false, ) + reset_containers!(containers) + # - Sample Function - # for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) # Define function to integrate @@ -48,6 +50,7 @@ function self_vortex_panel_integration!( containers; debug=false, ) + reset_containers!(containers) # - Sample Function - # for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) @@ -93,6 +96,7 @@ function nominal_source_panel_integration!( containers; debug=false, ) + reset_containers!(containers) # - Sample Function - # for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) @@ -127,6 +131,7 @@ function self_source_panel_integration!( containers; debug=false, ) + reset_containers!(containers) # - Sample Function - # for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) diff --git a/src/preprocess/velocities/induced_velocity_matrices.jl b/src/preprocess/velocities/induced_velocity_matrices.jl index d575b6fc..03dfb963 100644 --- a/src/preprocess/velocities/induced_velocity_matrices.jl +++ b/src/preprocess/velocities/induced_velocity_matrices.jl @@ -44,7 +44,7 @@ function induced_velocities_from_vortex_panels_on_points( influence_lengths, strengths, integration_options; - cache_vec=cache_vec, + integration_caches=integration_caches, ) return VEL @@ -170,7 +170,7 @@ function induced_velocities_from_source_panels_on_points( controlpoints, nodes, nodemap, - influence_lengthss, + influence_lengths, strengths, integration_options; integration_caches=integration_caches, diff --git a/src/preprocess/velocities/romberg_integrals.jl b/src/preprocess/velocities/romberg_integrals.jl index db5e63ab..ce141259 100644 --- a/src/preprocess/velocities/romberg_integrals.jl +++ b/src/preprocess/velocities/romberg_integrals.jl @@ -36,18 +36,17 @@ function extrapolate!(V, err, fh; power=2, atol=1e-6) for i in (numeval - 1):-1:1 f_i, h_i = fh[i + (firstindex(fh) - 1)] c = (h_i[] / h_prime[])^power - println((f_ip1 - f_i)) @. f_ip1 += (f_ip1 - f_i) / (c - 1) fh[i + (firstindex(fh) - 1)] = (f_ip1, h_i) err_prime = norm.(f_ip1 - f_i) minerr_prime = min.(minerr_prime, err_prime) if err_prime < err - V, err = f_ip1, err_prime + V[:], err = f_ip1, err_prime end end # - converged - # - all(err .<= atol) && break + # all(err .<= atol) && break end return V, err @@ -67,6 +66,7 @@ function nominal_vortex_panel_integration!( containers; debug=false, ) + reset_containers!(containers) # - Loop through number of subdivisions (start with 2) - # for ii in 1:(integration_options.max_subdivisions) @@ -82,7 +82,7 @@ function nominal_vortex_panel_integration!( t = (i - 0.5) / nint # get sample velocity - dt.nominal_vortex_induced_velocity_sample!( + nominal_vortex_induced_velocity_sample!( @view(containers.samples[ii])[1][1], t, node1, @@ -98,7 +98,7 @@ function nominal_vortex_panel_integration!( if ii > 1 # extrapolate once you can - dt.extrapolate!( + extrapolate!( V, @view(containers.sample_cache[1:4]), @view(containers.samples[1:ii]); @@ -123,10 +123,10 @@ function self_vortex_panel_integration!( node2, influence_length, controlpoint, - containers, + containers; debug=false, ) - + reset_containers!(containers) # - Loop through number of subdivisions (start with 2) - # for ii in 1:(integration_options.max_subdivisions) nint = 2^ii @@ -141,7 +141,7 @@ function self_vortex_panel_integration!( t = (i - 0.5) / nint # get sample velocity - dt.self_vortex_induced_velocity_sample!( + self_vortex_induced_velocity_sample!( @view(containers.samples[ii])[1][1], t, node1, @@ -165,7 +165,7 @@ function self_vortex_panel_integration!( if ii > 1 # extrapolate once you can - dt.extrapolate!( + extrapolate!( V, @view(containers.sample_cache[1:4]), @view(containers.samples[1:ii]); @@ -198,6 +198,7 @@ function nominal_source_panel_integration!( containers; debug=false, ) + reset_containers!(containers) # - Loop through number of subdivisions (start with 2) - # for ii in 1:(integration_options.max_subdivisions) @@ -213,7 +214,7 @@ function nominal_source_panel_integration!( t = (i - 0.5) / nint # get sample velocity - dt.nominal_source_induced_velocity_sample!( + nominal_source_induced_velocity_sample!( @view(containers.samples[ii])[1][1], t, node1, @@ -229,7 +230,7 @@ function nominal_source_panel_integration!( if ii > 1 # extrapolate once you can - dt.extrapolate!( + extrapolate!( V, @view(containers.sample_cache[1:4]), @view(containers.samples[1:ii]); @@ -257,7 +258,7 @@ function self_source_panel_integration!( containers; debug=false, ) - + reset_containers!(containers) # - Loop through number of subdivisions (start with 2) - # for ii in 1:(integration_options.max_subdivisions) nint = 2^ii @@ -272,7 +273,7 @@ function self_source_panel_integration!( t = (i - 0.5) / nint # get sample velocity - dt.self_source_induced_velocity_sample!( + self_source_induced_velocity_sample!( @view(containers.samples[ii])[1][1], t, node1, @@ -296,7 +297,7 @@ function self_source_panel_integration!( if ii > 1 # extrapolate once you can - dt.extrapolate!( + extrapolate!( V, @view(containers.sample_cache[1:4]), @view(containers.samples[1:ii]); diff --git a/src/utilities/types.jl b/src/utilities/types.jl index b3cc4038..dbcf90e9 100644 --- a/src/utilities/types.jl +++ b/src/utilities/types.jl @@ -156,19 +156,15 @@ end # QUADRATURE TYPES # #---------------------------------# -@kwdef struct Romberg{TI,TF,TP,TV} <: IntegrationMethod +@kwdef struct Romberg{TI,TF} <: IntegrationMethod max_subdivisions::TI = 10 atol::TF = 1e-6 end -function set_romberg_options(; max_subdivisions=10, atol=1e-6) - return Romberg(; max_subdivisions, atol) -end - @kwdef struct GaussKronrod{TI,TF} <: IntegrationMethod - order::TI = 5 - maxevals::TI = 1000 - atol::TF = 1e-12 + order::TI = 7 + maxevals::TI = 10^7 + atol::TF = 0 end struct GaussLegendre{TN,TW} <: IntegrationMethod diff --git a/src/utilities/utils.jl b/src/utilities/utils.jl index 998a7041..41155709 100644 --- a/src/utilities/utils.jl +++ b/src/utilities/utils.jl @@ -99,8 +99,16 @@ function reset_containers!(c) for p in propertynames(c) cp = getfield(c, p) if typeof(cp) <: AbstractArray - #do nothing if it's a string - (eltype(cp) == String) || (cp .= 0) + if eltype(cp) <: Tuple + for i in 1:length(cp[1]) + for j in 1:length(cp) + cp[j][i] .= 0.0 + end + end + else + #do nothing if it's a string + (eltype(cp) == String) || (cp .= 0) + end else reset_containers!(cp) end diff --git a/test/induced_velocities.jl b/test/induced_velocities.jl index a957098f..74ef235a 100644 --- a/test/induced_velocities.jl +++ b/test/induced_velocities.jl @@ -88,14 +88,45 @@ end # - Test 1 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities1.jl") + + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + node1 = p1 node2 = p2 influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; + ) + + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-6) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-6) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-6) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-6) + + # GaussLegendre + V = dt.nominal_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-6) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-6) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-6) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-6) + + # Romberg + V = dt.nominal_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; ) # Compare with DFDC integration values @@ -112,9 +143,10 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) # Compare with DFDC integration values @@ -123,6 +155,29 @@ end @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-9) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-9) + # GaussLegendre + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-9) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-9) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-9) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-9) + + # Romberg + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # - Test 3 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities3.jl") @@ -131,17 +186,40 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + + # GaussLegendre # Calculate Integral V = dt.nominal_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # - Test 4 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities4.jl") @@ -150,11 +228,33 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + + # Romberg + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @@ -164,6 +264,14 @@ end @testset "Single Vortex Panel Self-Induction Integration" begin + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + + # - Test 1 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities1.jl") @@ -172,16 +280,41 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) #note: radial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) @@ -194,20 +327,47 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) - # Compare with DFDC integration values # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) #note: radial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-3) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-3) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-3) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-3) + + + # - Test 3 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities3.jl") @@ -216,20 +376,45 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) - # Compare with DFDC integration values # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) - #note: radial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-4) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + + # Romberg + # Calculate Integral + V = dt.self_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-3) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-3) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-3) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-3) + # - Test 4 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities4.jl") @@ -238,27 +423,61 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) #note: radial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + end + #---------------------------------# # SOURCES # #---------------------------------# @testset "Nominal Single Source Panel Integration" begin + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + # - Test 1 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities1.jl") @@ -267,17 +486,39 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-6) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-6) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-6) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-6) + # GaussLegendre + # Calculate Integral + V = dt.nominal_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vsigmai[1], atol=1e-6) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-6) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-6) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-6) + # Romberg + # Calculate Integral + V = dt.nominal_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-4) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-4) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-4) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-4) + # - Test 2 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities2.jl") @@ -286,17 +527,41 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-9) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-9) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-9) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-9) + # GaussLegendre + # Calculate Integral + V = dt.nominal_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vsigmai[1], atol=1e-9) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-9) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-9) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-9) + # Romber + # Calculate Integral + V = dt.nominal_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-4) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-4) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-4) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-4) + + + # - Test 3 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities3.jl") @@ -305,17 +570,41 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.nominal_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.nominal_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-4) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-4) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-4) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-4) + + + # - Test 4 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities4.jl") @@ -324,20 +613,50 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.nominal_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + + # Romber + # Calculate Integral + V = dt.nominal_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + end @testset "Single Source Panel Self-Induction Integration" begin + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + # - Test 1 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities1.jl") @@ -346,19 +665,44 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values #note: axial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + + # - Test 2 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities2.jl") @@ -367,19 +711,44 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) - # Compare with DFDC integration values #note: axial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-4) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-4) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-4) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-4) + + + # - Test 3 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities3.jl") @@ -388,19 +757,42 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) - # Compare with DFDC integration values #note: axial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-3) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-3) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-3) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-3) + # - Test 4 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities4.jl") @@ -409,21 +801,45 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values #note: axial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) end end + ###################################################################### # # # Multi-PANEL INDUCED VELOCITIES # @@ -432,6 +848,14 @@ end @testset "Multi-Panel, Multi-Target Induced Velocity Matrices" begin @testset "Multiple Vortex Panel Induced Velocities on Multiple Points" begin + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + + # define control points controlpoints = [0.0 1.0; 1.0 1.0]' @@ -447,26 +871,52 @@ end # use unit strengths strengths = [1.0 1.0; 1.0 1.0] + # - GaussKronrod - # # [cp, n, x/r] VEL = dt.induced_velocities_from_vortex_panels_on_points( - controlpoints, nodes, nodemap, influence_lengths, strengths + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=gk_integration_options, singular=gk_integration_options), ) # [vz1 vr1; vz2 vr2] vn12cp1 = dt.self_vortex_panel_integration( - nodes[:, 1], nodes[:, 2], influence_lengths[1], controlpoints[:, 1], zeros(20) + gk_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + gk_cache, ) vn12cp2 = dt.nominal_vortex_panel_integration( - nodes[:, 1], nodes[:, 2], influence_lengths[1], controlpoints[:, 2], zeros(20) + gk_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + gk_cache, ) vn23cp1 = dt.nominal_vortex_panel_integration( - nodes[:, 2], nodes[:, 3], influence_lengths[1], controlpoints[:, 1], zeros(20) + gk_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + gk_cache, ) vn23cp2 = dt.nominal_vortex_panel_integration( - nodes[:, 2], nodes[:, 3], influence_lengths[1], controlpoints[:, 2], zeros(20) + gk_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + gk_cache, ) @test all(VEL[1, 1, :] .== vn12cp1[1, :]) @@ -476,10 +926,131 @@ end @test all(VEL[2, 1, :] .== vn12cp2[1, :]) @test all(VEL[2, 2, :] .== vn12cp2[2, :] .+ vn23cp2[1, :]) @test all(VEL[2, 3, :] .== vn23cp2[2, :]) + + # - GaussLegendre - # + # [cp, n, x/r] + VEL = dt.induced_velocities_from_vortex_panels_on_points( + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=gl_integration_options, singular=gl_integration_options), + ) + + # [vz1 vr1; vz2 vr2] + vn12cp1 = dt.self_vortex_panel_integration( + gl_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + gl_cache, + ) + + vn12cp2 = dt.nominal_vortex_panel_integration( + gl_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + gl_cache, + ) + + vn23cp1 = dt.nominal_vortex_panel_integration( + gl_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + gl_cache, + ) + + vn23cp2 = dt.nominal_vortex_panel_integration( + gl_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + gl_cache, + ) + + @test all(VEL[1, 1, :] .== vn12cp1[1, :]) + #note that node connected to 2 panels has influence contributions from both panels + @test all(VEL[1, 2, :] .== vn12cp1[2, :] .+ vn23cp1[1, :]) + @test all(VEL[1, 3, :] .== vn23cp1[2, :]) + @test all(VEL[2, 1, :] .== vn12cp2[1, :]) + @test all(VEL[2, 2, :] .== vn12cp2[2, :] .+ vn23cp2[1, :]) + @test all(VEL[2, 3, :] .== vn23cp2[2, :]) + + # - Romberg - # + # [cp, n, x/r] + VEL = dt.induced_velocities_from_vortex_panels_on_points( + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=r_integration_options, singular=r_integration_options), + ) + + # [vz1 vr1; vz2 vr2] + vn12cp1 = dt.self_vortex_panel_integration( + r_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + r_cache, + ) + + vn12cp2 = dt.nominal_vortex_panel_integration( + r_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + r_cache, + ) + + vn23cp1 = dt.nominal_vortex_panel_integration( + r_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + r_cache, + ) + + vn23cp2 = dt.nominal_vortex_panel_integration( + r_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + r_cache, + ) + + @test all(VEL[1, 1, :] .== vn12cp1[1, :]) + #note that node connected to 2 panels has influence contributions from both panels + @test all(VEL[1, 2, :] .== vn12cp1[2, :] .+ vn23cp1[1, :]) + @test all(VEL[1, 3, :] .== vn23cp1[2, :]) + #TODO: why is the answer different when directly vs indirectly calling integration? + @test all(isapprox.(VEL[2, 1, :] , vn12cp2[1, :], atol=1e-4)) + @test all(isapprox.(VEL[2, 2, :] , vn12cp2[2, :] .+ vn23cp2[1, :], atol=1e-4)) + @test all(VEL[2, 3, :] .== vn23cp2[2, :]) end @testset "Multiple Source Panel Induced Velocities on Multiple Points" begin + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + + # define control points controlpoints = [0.0 1.0; 1.0 1.0]' @@ -495,26 +1066,107 @@ end # use unit strengths strengths = [1.0 1.0; 1.0 1.0] + # GaussKronrod + # [cp, n, x/r] + VEL = dt.induced_velocities_from_source_panels_on_points( + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=gk_integration_options, singular=gk_integration_options), + ) + + # [vz1 vr1; vz2 vr2] + vn12cp1 = dt.self_source_panel_integration( + gk_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + gk_cache, + ) + + vn12cp2 = dt.nominal_source_panel_integration( + gk_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + gk_cache, + ) + + vn23cp1 = dt.nominal_source_panel_integration( + gk_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + gk_cache, + ) + + vn23cp2 = dt.nominal_source_panel_integration( + gk_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + gk_cache, + ) + + @test all(VEL[1, 1, :] .== vn12cp1[1, :]) + @test all(VEL[1, 2, :] .== vn12cp1[2, :] .+ vn23cp1[1, :]) + @test all(VEL[1, 3, :] .== vn23cp1[2, :]) + @test all(VEL[2, 1, :] .== vn12cp2[1, :]) + @test all(VEL[2, 2, :] .== vn12cp2[2, :] .+ vn23cp2[1, :]) + @test all(VEL[2, 3, :] .== vn23cp2[2, :]) + + # GaussLegendre # [cp, n, x/r] VEL = dt.induced_velocities_from_source_panels_on_points( - controlpoints, nodes, nodemap, influence_lengths, strengths + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=gl_integration_options, singular=gl_integration_options), ) # [vz1 vr1; vz2 vr2] vn12cp1 = dt.self_source_panel_integration( - nodes[:, 1], nodes[:, 2], influence_lengths[1], controlpoints[:, 1], zeros(20) + gl_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + gl_cache, ) vn12cp2 = dt.nominal_source_panel_integration( - nodes[:, 1], nodes[:, 2], influence_lengths[1], controlpoints[:, 2], zeros(20) + gl_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + gl_cache, ) vn23cp1 = dt.nominal_source_panel_integration( - nodes[:, 2], nodes[:, 3], influence_lengths[1], controlpoints[:, 1], zeros(20) + gl_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + gl_cache, ) vn23cp2 = dt.nominal_source_panel_integration( - nodes[:, 2], nodes[:, 3], influence_lengths[1], controlpoints[:, 2], zeros(20) + gl_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + gl_cache, ) @test all(VEL[1, 1, :] .== vn12cp1[1, :]) @@ -523,5 +1175,61 @@ end @test all(VEL[2, 1, :] .== vn12cp2[1, :]) @test all(VEL[2, 2, :] .== vn12cp2[2, :] .+ vn23cp2[1, :]) @test all(VEL[2, 3, :] .== vn23cp2[2, :]) + + # Romberg + # [cp, n, x/r] + VEL = dt.induced_velocities_from_source_panels_on_points( + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=r_integration_options, singular=r_integration_options), + ) + + # [vz1 vr1; vz2 vr2] + vn12cp1 = dt.self_source_panel_integration( + r_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + r_cache, + ) + + vn12cp2 = dt.nominal_source_panel_integration( + r_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + r_cache, + ) + + vn23cp1 = dt.nominal_source_panel_integration( + r_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + r_cache, + ) + + vn23cp2 = dt.nominal_source_panel_integration( + r_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + r_cache, + ) + + @test all(VEL[1, 1, :] .== vn12cp1[1, :]) + @test all(VEL[1, 2, :] .== vn12cp1[2, :] .+ vn23cp1[1, :]) + @test all(VEL[1, 3, :] .== vn23cp1[2, :]) + # TODO: why are these specific two different (but not totally wrong)?? + @test all(isapprox.(VEL[2, 1, :] , vn12cp2[1, :], atol=1e-4)) + @test all(isapprox.(VEL[2, 2, :] , vn12cp2[2, :] .+ vn23cp2[1, :],atol=1e-3)) + @test all(VEL[2, 3, :] .== vn23cp2[2, :]) end end diff --git a/test/linear_system_assembly.jl b/test/linear_system_assembly.jl index a175565d..42b52a5d 100644 --- a/test/linear_system_assembly.jl +++ b/test/linear_system_assembly.jl @@ -1,5 +1,6 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") @testset "Linear System Assembly" begin + integration_options = (; nominal=dt.GaussLegendre(), singular=dt.GaussLegendre()) # define coordinates x1 = [1.0; 0.5; 0.0; 0.5; 1.0] @@ -25,6 +26,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.node, panels.nodemap, panels.influence_length, + integration_options, ) dt.add_te_gap_aic!( @@ -36,6 +38,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.tendotn, panels.tencrossn, panels.teadjnodeidxs, + integration_options ) AICpcp = dt.vortex_aic_boundary_on_field( @@ -44,6 +47,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.node, panels.nodemap, panels.influence_length, + integration_options ) dt.add_te_gap_aic!( @@ -55,6 +59,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.tendotn, panels.tencrossn, panels.teadjnodeidxs, + integration_options ) Vinf = 1.0 #magnitude doesn't matter yet. @@ -129,6 +134,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.node, panels.nodemap, panels.influence_length, + integration_options ) AICpcp = dt.vortex_aic_boundary_on_field( @@ -137,6 +143,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.node, panels.nodemap, panels.influence_length, + integration_options ) Vinf = 1.0 #magnitude doesn't matter yet. diff --git a/test/runtests.jl b/test/runtests.jl index 5e1aac8f..2f724694 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,16 +23,13 @@ include("test_utils.jl") println("Running Tests...") -#---------------------------------# -# New Tests # -#---------------------------------# - # - Caching Tests - # include("iotests.jl") # - pre-process related tests - # include("afcorrections.jl") include("panel_generation_tests.jl") +include("induced_velocities.jl") include("influence_coefficients.jl") include("linear_system_assembly.jl") include("pre_processing_tests.jl") From ccee801d234f0ae483669ad5cfe25917207d8274 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Thu, 4 Apr 2024 10:43:47 -0600 Subject: [PATCH 12/12] update convenience function --- src/analysis/analyses.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/analysis/analyses.jl b/src/analysis/analyses.jl index 53187c64..396dd90a 100644 --- a/src/analysis/analyses.jl +++ b/src/analysis/analyses.jl @@ -91,7 +91,7 @@ function analyze( solve_parameter_tuple.linsys, solve_parameter_tuple.wakeK, propulsor; - grid_solver_options=options.wake_solver_options, + grid_solver_options=options.grid_solver_options, integration_options=options.integration_options, autoshiftduct=options.autoshiftduct, itcpshift=options.itcpshift, @@ -106,11 +106,11 @@ function analyze( #= NOTE: If the linear system or wake did not converge, there is likely a serious problem that would lead to an error in the solve, so we will exit here with a fail flag for an optimizer or user =# - if iszero(lu_decomp_flag) || !options.wake_solver_options.converged[1] + if iszero(lu_decomp_flag) || !options.grid_solver_options.converged[1] if !options.silence_warnings if iszero(lu_decomp_flag) @warn "Exiting. LU decomposition of the LHS matrix for the linear system failed. Please check your body geometry and ensure that there will be no panels lying directly atop eachother or other similar problematic geometry." - elseif !options.wake_solver_options.converged[1] + elseif !options.grid_solver_options.converged[1] @warn "Exiting. Wake elliptic grid solve did not converge. Consider a looser convergence tolerance if the geometry looks good." end end