diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index 6ed0063c..989e5ec8 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 @@ -86,12 +90,16 @@ 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/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") ##### ----- PROCESS ----- ##### - # Solve and Residual Functions include("process/solve.jl") include("process/residuals/CSORresidual.jl") diff --git a/src/analysis/analyses.jl b/src/analysis/analyses.jl index 2923873e..396dd90a 100644 --- a/src/analysis/analyses.jl +++ b/src/analysis/analyses.jl @@ -84,14 +84,15 @@ 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, 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, axistol=options.axistol, @@ -105,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 diff --git a/src/preprocess/preprocess.jl b/src/preprocess/preprocess.jl index 858c5bd3..924c692e 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, @@ -1539,14 +1580,15 @@ function precompute_parameters!( operating_point, reference_parameters, problem_dimensions=nothing; - grid_solver_options=options.wake_solver_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 - # @@ -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..53f90079 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(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 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_caches, ) 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 new file mode 100644 index 00000000..d547ff70 --- /dev/null +++ b/src/preprocess/velocities/gausskronrod_integrals.jl @@ -0,0 +1,171 @@ +#---------------------------------# +# 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, + sample_cache; + debug=false, +) + + # Define function to integrate + function fsample(t) + return nominal_vortex_induced_velocity_sample( + t, node1, node2, influence_length, controlpoint, sample_cache + ) + 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( + integration_options::GaussKronrod, + node1, + node2, + influence_length, + controlpoint, + sample_cache; + debug=false, +) + + # Define function to integrate + function fsample(t) + return self_vortex_induced_velocity_sample( + t, node1, node2, influence_length, controlpoint, sample_cache; + ) + end + + V, err = quadgk( + fsample, + 0.0, + 0.5, + 1.0; + order=integration_options.order, + maxevals=integration_options.maxevals, + atol=integration_options.atol, + ) + + sample_cache[1], sample_cache[2] = analytically_integrated_vortex_influence( + controlpoint[2], influence_length + ) + + V .*= influence_length + V[1:2] .+= sample_cache[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( + integration_options::GaussKronrod, + node1, + node2, + influence_length, + controlpoint, + sample_cache; + debug=false, +) + + # Define function to integrate + function fsample(t) + return nominal_source_induced_velocity_sample( + t, node1, node2, influence_length, controlpoint, sample_cache; + ) + 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( + integration_options::GaussKronrod, + node1, + node2, + influence_length, + controlpoint, + sample_cache; + debug=false, +) + + # Define function to integrate + function fsample(t) + return self_source_induced_velocity_sample( + t, node1, node2, influence_length, controlpoint, sample_cache; + ) + end + + V, err = quadgk( + fsample, + 0.0, + 0.5, + 1.0; + order=integration_options.order, + maxevals=integration_options.maxevals, + atol=integration_options.atol, + ) + + sample_cache[1], sample_cache[2] = analytically_integrated_source_influence( + controlpoint[2], influence_length + ) + + V .*= influence_length + V[3:4] .+= 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/preprocess/velocities/gausslegendre_integrals.jl b/src/preprocess/velocities/gausslegendre_integrals.jl new file mode 100644 index 00000000..1c1ff6a0 --- /dev/null +++ b/src/preprocess/velocities/gausslegendre_integrals.jl @@ -0,0 +1,161 @@ +#---------------------------------# +# 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, +) + reset_containers!(containers) + + # - Sample Function - # + for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) + # Define function to integrate + 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 + +""" +`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, +) + reset_containers!(containers) + + # - 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 # +#---------------------------------# + +""" +`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, +) + reset_containers!(containers) + + # - Sample Function - # + for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) + # Define function to integrate + 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 + +""" +`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, +) + reset_containers!(containers) + + # - 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/preprocess/velocities/induced_velocity_matrices.jl b/src/preprocess/velocities/induced_velocity_matrices.jl index 949c3d17..03dfb963 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,8 +42,9 @@ function induced_velocities_from_vortex_panels_on_points( nodes, nodemap, influence_lengths, - strengths; - cache_vec=cache_vec, + strengths, + integration_options; + integration_caches=integration_caches, ) return VEL @@ -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 @@ -133,8 +171,9 @@ function induced_velocities_from_source_panels_on_points( nodes, nodemap, influence_lengths, - strengths; - cache_vec=cache_vec, + 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/integrals.jl b/src/preprocess/velocities/integrands.jl similarity index 57% rename from src/preprocess/velocities/integrals.jl rename to src/preprocess/velocities/integrands.jl index d7f75b33..4511a836 100644 --- a/src/preprocess/velocities/integrals.jl +++ b/src/preprocess/velocities/integrands.jl @@ -1,12 +1,54 @@ ###################################################################### # # -# Linear Vortex Panel Integration # +# VORTEX # # # ###################################################################### #---------------------------------# # 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] = 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)) + + # 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 """ @@ -18,8 +60,8 @@ function nominal_vortex_induced_velocity_sample( ) # Transform from (0,1) to actual position on panel - cache_vec[1] = FLOWMath.linear(nondimrange, (node1[1], node2[1]), t) # z coordinate - cache_vec[2] = FLOWMath.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)) @@ -58,43 +100,6 @@ function nominal_vortex_induced_velocity_sample( # 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}] -""" -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 # #---------------------------------# @@ -147,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) @@ -156,17 +169,16 @@ 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( - t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] +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] = FLOWMath.linear(nondimrange, (node1[1], node2[1]), t) - cache_vec[2] = FLOWMath.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)) @@ -195,59 +207,63 @@ function self_vortex_induced_velocity_sample( 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[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 """ -`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}] +# Arguments: +- `t::Float` : sample point in range (0,1) selected by quadrature. """ -function self_vortex_panel_integration( - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=[0.0; 1.0], - debug=false, +function self_vortex_induced_velocity_sample( + t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] ) - # Define function to integrate - function fsample(t) - return self_vortex_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, - ) - end + # 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) - V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) + # get relative geometry + calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) - cache_vec[1], cache_vec[2] = analytically_integrated_vortex_influence( - controlpoint[2], influence_length - ) + # 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 - V .*= influence_length - V[1:2] .+= cache_vec[1] / 2.0 + # 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) + ) - if debug - return reshape(V, (2, 2)), err - else - return reshape(V, (2, 2)) - end + #= + 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 ###################################################################### # # -# Linear Source Panel Integration # +# SOURCE # # # ###################################################################### @@ -256,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. """ @@ -265,8 +320,8 @@ function nominal_source_induced_velocity_sample( ) # Transform from (0,1) to actual position on panel - cache_vec[1] = FLOWMath.linear(nondimrange, (node1[1], node2[1]), t) - cache_vec[2] = FLOWMath.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)) @@ -300,43 +355,6 @@ function nominal_source_induced_velocity_sample( ) 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( - 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 # #---------------------------------# @@ -375,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) @@ -384,17 +410,16 @@ 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( - t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] +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] = FLOWMath.linear(nondimrange, (node1[1], node2[1]), t) - cache_vec[2] = FLOWMath.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)) @@ -422,52 +447,56 @@ function self_source_induced_velocity_sample( 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 - ] + 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 """ -`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}] +# Arguments: +- `t::Float` : sample point in range (0,1) selected by quadrature. """ -function self_source_panel_integration( - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=[0.0; 1.0], - debug=false, +function self_source_induced_velocity_sample( + t, node1, node2, influence_length, controlpoint, cache_vec; nondimrange=[0.0; 1.0] ) - # Define function to integrate - function fsample(t) - return self_source_induced_velocity_sample( - t, - node1, - node2, - influence_length, - controlpoint, - cache_vec; - nondimrange=nondimrange, - ) - end + # 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) - V, err = quadgk(fsample, 0.0, 0.5, 1.0; order=3, maxevals=1e2, atol=1e-6) + # get relative geometry + calculate_xrm!(view(cache_vec, 3:6), controlpoint, view(cache_vec, 1:2)) - cache_vec[1], cache_vec[2] = analytically_integrated_source_influence( - controlpoint[2], influence_length + # 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) ) - V .*= influence_length - V[3:4] .+= cache_vec[2] / 2.0 + # 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) + ) - if debug - return reshape(V, (2, 2)), err - else - return reshape(V, (2, 2)) - end + #= + 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/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 new file mode 100644 index 00000000..ce141259 --- /dev/null +++ b/src/preprocess/velocities/romberg_integrals.jl @@ -0,0 +1,318 @@ +""" + 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!(V, err, fh; power=2, atol=1e-6) + + # - rename for convenience - # + (f0, h) = first(fh) + + # - initialize - # + err .= Inf + numeval = 1 + maxeval = size(fh, 1) + + # - Loop - # + while numeval < maxeval + + # - increment - # + numeval += 1 + + # - update - # + (f_prime, h_prime) = fh[numeval + (firstindex(fh) - 1)] + h = h_prime + 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) + 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 + end + end + + # - converged - # + # all(err .<= atol) && break + end + + return V, err +end + +#---------------------------------# +# VORTEX # +#---------------------------------# + +function nominal_vortex_panel_integration!( + integration_options::Romberg, + V, + node1, + node2, + influence_length, + controlpoint, + 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 + + # - 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 + 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 + 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 self_vortex_panel_integration!( + integration_options::Romberg, + V, + node1, + node2, + influence_length, + controlpoint, + 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 + + # - 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 + 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 + 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 V, err + else + return V + end +end + +#---------------------------------# +# SOURCE # +#---------------------------------# + +function nominal_source_panel_integration!( + integration_options::Romberg, + V, + node1, + node2, + influence_length, + controlpoint, + 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 + + # - 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 + 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 + 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 self_source_panel_integration!( + integration_options::Romberg, + V, + node1, + node2, + influence_length, + controlpoint, + 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 + + # - 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 + 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 + 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 V, err + else + return V + end +end + 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 7eb70fb0..dbcf90e9 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,41 @@ function verify_input(propulsor) # TODO: go find all the various asserts and put them here end +#---------------------------------# +# QUADRATURE TYPES # +#---------------------------------# + +@kwdef struct Romberg{TI,TF} <: IntegrationMethod + max_subdivisions::TI = 10 + atol::TF = 1e-6 +end + +@kwdef struct GaussKronrod{TI,TF} <: IntegrationMethod + order::TI = 7 + maxevals::TI = 10^7 + atol::TF = 0 +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 +264,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 +278,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" 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")