diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index 82db5319..989e5ec8 100644 --- a/src/DuctAPE.jl +++ b/src/DuctAPE.jl @@ -41,8 +41,6 @@ using ForwardDiff # used for jacobian for newton solver using FLOWMath # used for various items, mostly interpolation using Printf # used when verbose option is selected -using Primes # used in Romberg integration settings - #---------------------------------# # EXPORTS # #---------------------------------# diff --git a/src/preprocess/preprocess.jl b/src/preprocess/preprocess.jl index 5e7d86fa..924c692e 100644 --- a/src/preprocess/preprocess.jl +++ b/src/preprocess/preprocess.jl @@ -1580,15 +1580,15 @@ function precompute_parameters!( operating_point, reference_parameters, problem_dimensions=nothing; - grid_solver_options=options.wake_solver_options, - integration_options=options.integration_options, - autoshiftduct=options.autoshiftduct, - itcpshift=options.itcpshift, - axistol=options.axistol, - tegaptol=options.tegaptol, - finterp=options.finterp, - silence_warnings=options.silence_warnings, - verbose=options.verbose, + grid_solver_options=GridSolverOptions(), + integration_options=IntegrationOptions(), + autoshiftduct=true, + itcpshift=0.05, + axistol=1e-15, + tegaptol=1e1 * eps(), + finterp=fm.akima, + silence_warnings=true, + verbose=false, ) # - Reset Caches - # diff --git a/src/preprocess/velocities/body_aic.jl b/src/preprocess/velocities/body_aic.jl index 73091a11..53f90079 100644 --- a/src/preprocess/velocities/body_aic.jl +++ b/src/preprocess/velocities/body_aic.jl @@ -72,10 +72,10 @@ function vortex_aic_boundary_on_boundary!( if isnothing(integration_caches) # integration_cache = zeros(eltype(controlpoint), 20) nominal_integration_cache = allocate_integration_containers( - integration_options.nominal, eltype(VEL) + integration_options.nominal, eltype(AICn) ) singular_integration_cache = allocate_integration_containers( - integration_options.singular, eltype(VEL) + integration_options.singular, eltype(AICn) ) else nominal_integration_cache = integration_caches.nominal @@ -164,7 +164,7 @@ function vortex_aic_boundary_on_field( nodemap, influence_length, integration_options; - integration_caches=integration_cache, + integration_caches=integration_caches, ) return AICn diff --git a/src/preprocess/velocities/gausslegendre_integrals.jl b/src/preprocess/velocities/gausslegendre_integrals.jl index 546e4141..1c1ff6a0 100644 --- a/src/preprocess/velocities/gausslegendre_integrals.jl +++ b/src/preprocess/velocities/gausslegendre_integrals.jl @@ -15,6 +15,8 @@ function nominal_vortex_panel_integration!( containers; debug=false, ) + reset_containers!(containers) + # - Sample Function - # for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) # Define function to integrate @@ -48,6 +50,7 @@ function self_vortex_panel_integration!( containers; debug=false, ) + reset_containers!(containers) # - Sample Function - # for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) @@ -93,6 +96,7 @@ function nominal_source_panel_integration!( containers; debug=false, ) + reset_containers!(containers) # - Sample Function - # for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) @@ -127,6 +131,7 @@ function self_source_panel_integration!( containers; debug=false, ) + reset_containers!(containers) # - Sample Function - # for (s, t) in zip(eachrow(containers.samples), integration_options.sample_points) diff --git a/src/preprocess/velocities/induced_velocity_matrices.jl b/src/preprocess/velocities/induced_velocity_matrices.jl index d575b6fc..03dfb963 100644 --- a/src/preprocess/velocities/induced_velocity_matrices.jl +++ b/src/preprocess/velocities/induced_velocity_matrices.jl @@ -44,7 +44,7 @@ function induced_velocities_from_vortex_panels_on_points( influence_lengths, strengths, integration_options; - cache_vec=cache_vec, + integration_caches=integration_caches, ) return VEL @@ -170,7 +170,7 @@ function induced_velocities_from_source_panels_on_points( controlpoints, nodes, nodemap, - influence_lengthss, + influence_lengths, strengths, integration_options; integration_caches=integration_caches, diff --git a/src/preprocess/velocities/romberg_integrals.jl b/src/preprocess/velocities/romberg_integrals.jl index db5e63ab..ce141259 100644 --- a/src/preprocess/velocities/romberg_integrals.jl +++ b/src/preprocess/velocities/romberg_integrals.jl @@ -36,18 +36,17 @@ function extrapolate!(V, err, fh; power=2, atol=1e-6) for i in (numeval - 1):-1:1 f_i, h_i = fh[i + (firstindex(fh) - 1)] c = (h_i[] / h_prime[])^power - println((f_ip1 - f_i)) @. f_ip1 += (f_ip1 - f_i) / (c - 1) fh[i + (firstindex(fh) - 1)] = (f_ip1, h_i) err_prime = norm.(f_ip1 - f_i) minerr_prime = min.(minerr_prime, err_prime) if err_prime < err - V, err = f_ip1, err_prime + V[:], err = f_ip1, err_prime end end # - converged - # - all(err .<= atol) && break + # all(err .<= atol) && break end return V, err @@ -67,6 +66,7 @@ function nominal_vortex_panel_integration!( containers; debug=false, ) + reset_containers!(containers) # - Loop through number of subdivisions (start with 2) - # for ii in 1:(integration_options.max_subdivisions) @@ -82,7 +82,7 @@ function nominal_vortex_panel_integration!( t = (i - 0.5) / nint # get sample velocity - dt.nominal_vortex_induced_velocity_sample!( + nominal_vortex_induced_velocity_sample!( @view(containers.samples[ii])[1][1], t, node1, @@ -98,7 +98,7 @@ function nominal_vortex_panel_integration!( if ii > 1 # extrapolate once you can - dt.extrapolate!( + extrapolate!( V, @view(containers.sample_cache[1:4]), @view(containers.samples[1:ii]); @@ -123,10 +123,10 @@ function self_vortex_panel_integration!( node2, influence_length, controlpoint, - containers, + containers; debug=false, ) - + reset_containers!(containers) # - Loop through number of subdivisions (start with 2) - # for ii in 1:(integration_options.max_subdivisions) nint = 2^ii @@ -141,7 +141,7 @@ function self_vortex_panel_integration!( t = (i - 0.5) / nint # get sample velocity - dt.self_vortex_induced_velocity_sample!( + self_vortex_induced_velocity_sample!( @view(containers.samples[ii])[1][1], t, node1, @@ -165,7 +165,7 @@ function self_vortex_panel_integration!( if ii > 1 # extrapolate once you can - dt.extrapolate!( + extrapolate!( V, @view(containers.sample_cache[1:4]), @view(containers.samples[1:ii]); @@ -198,6 +198,7 @@ function nominal_source_panel_integration!( containers; debug=false, ) + reset_containers!(containers) # - Loop through number of subdivisions (start with 2) - # for ii in 1:(integration_options.max_subdivisions) @@ -213,7 +214,7 @@ function nominal_source_panel_integration!( t = (i - 0.5) / nint # get sample velocity - dt.nominal_source_induced_velocity_sample!( + nominal_source_induced_velocity_sample!( @view(containers.samples[ii])[1][1], t, node1, @@ -229,7 +230,7 @@ function nominal_source_panel_integration!( if ii > 1 # extrapolate once you can - dt.extrapolate!( + extrapolate!( V, @view(containers.sample_cache[1:4]), @view(containers.samples[1:ii]); @@ -257,7 +258,7 @@ function self_source_panel_integration!( containers; debug=false, ) - + reset_containers!(containers) # - Loop through number of subdivisions (start with 2) - # for ii in 1:(integration_options.max_subdivisions) nint = 2^ii @@ -272,7 +273,7 @@ function self_source_panel_integration!( t = (i - 0.5) / nint # get sample velocity - dt.self_source_induced_velocity_sample!( + self_source_induced_velocity_sample!( @view(containers.samples[ii])[1][1], t, node1, @@ -296,7 +297,7 @@ function self_source_panel_integration!( if ii > 1 # extrapolate once you can - dt.extrapolate!( + extrapolate!( V, @view(containers.sample_cache[1:4]), @view(containers.samples[1:ii]); diff --git a/src/utilities/types.jl b/src/utilities/types.jl index b3cc4038..dbcf90e9 100644 --- a/src/utilities/types.jl +++ b/src/utilities/types.jl @@ -156,19 +156,15 @@ end # QUADRATURE TYPES # #---------------------------------# -@kwdef struct Romberg{TI,TF,TP,TV} <: IntegrationMethod +@kwdef struct Romberg{TI,TF} <: IntegrationMethod max_subdivisions::TI = 10 atol::TF = 1e-6 end -function set_romberg_options(; max_subdivisions=10, atol=1e-6) - return Romberg(; max_subdivisions, atol) -end - @kwdef struct GaussKronrod{TI,TF} <: IntegrationMethod - order::TI = 5 - maxevals::TI = 1000 - atol::TF = 1e-12 + order::TI = 7 + maxevals::TI = 10^7 + atol::TF = 0 end struct GaussLegendre{TN,TW} <: IntegrationMethod diff --git a/src/utilities/utils.jl b/src/utilities/utils.jl index 998a7041..41155709 100644 --- a/src/utilities/utils.jl +++ b/src/utilities/utils.jl @@ -99,8 +99,16 @@ function reset_containers!(c) for p in propertynames(c) cp = getfield(c, p) if typeof(cp) <: AbstractArray - #do nothing if it's a string - (eltype(cp) == String) || (cp .= 0) + if eltype(cp) <: Tuple + for i in 1:length(cp[1]) + for j in 1:length(cp) + cp[j][i] .= 0.0 + end + end + else + #do nothing if it's a string + (eltype(cp) == String) || (cp .= 0) + end else reset_containers!(cp) end diff --git a/test/induced_velocities.jl b/test/induced_velocities.jl index a957098f..74ef235a 100644 --- a/test/induced_velocities.jl +++ b/test/induced_velocities.jl @@ -88,14 +88,45 @@ end # - Test 1 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities1.jl") + + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + node1 = p1 node2 = p2 influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; + ) + + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-6) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-6) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-6) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-6) + + # GaussLegendre + V = dt.nominal_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-6) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-6) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-6) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-6) + + # Romberg + V = dt.nominal_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; ) # Compare with DFDC integration values @@ -112,9 +143,10 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) # Compare with DFDC integration values @@ -123,6 +155,29 @@ end @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-9) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-9) + # GaussLegendre + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-9) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-9) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-9) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-9) + + # Romberg + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # - Test 3 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities3.jl") @@ -131,17 +186,40 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + + # GaussLegendre # Calculate Integral V = dt.nominal_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # - Test 4 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities4.jl") @@ -150,11 +228,33 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + + # Romberg + # Calculate Integral + V = dt.nominal_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @@ -164,6 +264,14 @@ end @testset "Single Vortex Panel Self-Induction Integration" begin + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + + # - Test 1 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities1.jl") @@ -172,16 +280,41 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) #note: radial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) @@ -194,20 +327,47 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) - # Compare with DFDC integration values # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) #note: radial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-3) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-3) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-3) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-3) + + + # - Test 3 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities3.jl") @@ -216,20 +376,45 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) - # Compare with DFDC integration values # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) - #note: radial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-4) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + + # Romberg + # Calculate Integral + V = dt.self_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-3) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-3) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-3) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-3) + # - Test 4 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities4.jl") @@ -238,27 +423,61 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_vortex_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_vortex_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) + #note: radial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) + @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_vortex_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + # note: axial tests fail if I use log(8Δs/r) in the analytic term rather than the 16 DFDC has instead of the 8 + @test isapprox(V[1, 1], Vgammai[1], atol=1e-5) + @test isapprox(V[2, 1], Vgammaip1[1], atol=1e-5) #note: radial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 2], Vgammai[2], atol=1e-5) @test isapprox(V[2, 2], Vgammaip1[2], atol=1e-5) + end + #---------------------------------# # SOURCES # #---------------------------------# @testset "Nominal Single Source Panel Integration" begin + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + # - Test 1 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities1.jl") @@ -267,17 +486,39 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-6) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-6) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-6) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-6) + # GaussLegendre + # Calculate Integral + V = dt.nominal_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vsigmai[1], atol=1e-6) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-6) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-6) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-6) + # Romberg + # Calculate Integral + V = dt.nominal_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-4) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-4) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-4) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-4) + # - Test 2 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities2.jl") @@ -286,17 +527,41 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-9) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-9) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-9) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-9) + # GaussLegendre + # Calculate Integral + V = dt.nominal_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vsigmai[1], atol=1e-9) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-9) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-9) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-9) + # Romber + # Calculate Integral + V = dt.nominal_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-4) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-4) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-4) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-4) + + + # - Test 3 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities3.jl") @@ -305,17 +570,41 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.nominal_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.nominal_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-4) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-4) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-4) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-4) + + + # - Test 4 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/nominal_velocities4.jl") @@ -324,20 +613,50 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = pf + # GaussKronrod # Calculate Integral V = dt.nominal_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.nominal_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + + # Romber + # Calculate Integral + V = dt.nominal_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + end @testset "Single Source Panel Self-Induction Integration" begin + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + # - Test 1 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities1.jl") @@ -346,19 +665,44 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values #note: axial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + + # - Test 2 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities2.jl") @@ -367,19 +711,44 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) - # Compare with DFDC integration values #note: axial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-4) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-4) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-4) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-4) + + + # - Test 3 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities3.jl") @@ -388,19 +757,42 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) - # Compare with DFDC integration values #note: axial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-3) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-3) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-3) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-3) + # - Test 4 - # # Load in comparision values from DFDC extraction include("./data/single_linear_panel_integration/self_velocities4.jl") @@ -409,21 +801,45 @@ end influence_length = sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2) controlpoint = ps + # GaussKronrod # Calculate Integral V = dt.self_source_panel_integration( - node1, node2, influence_length, controlpoint, zeros(20); nondimrange=[0.0; 1.0] + gk_integration_options, node1, node2, influence_length, controlpoint, gk_cache; ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # GaussLegendre + # Calculate Integral + V = dt.self_source_panel_integration( + gl_integration_options, node1, node2, influence_length, controlpoint, gl_cache; + ) # Compare with DFDC integration values #note: axial terms have zero for the analytic addition, so they work fine @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) + @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) + @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) + # Romberg + # Calculate Integral + V = dt.self_source_panel_integration( + r_integration_options, node1, node2, influence_length, controlpoint, r_cache; + ) + # Compare with DFDC integration values + #note: axial terms have zero for the analytic addition, so they work fine + @test isapprox(V[1, 1], Vsigmai[1], atol=1e-5) + @test isapprox(V[2, 1], Vsigmaip1[1], atol=1e-5) @test isapprox(V[1, 2], Vsigmai[2], atol=1e-5) @test isapprox(V[2, 2], Vsigmaip1[2], atol=1e-5) end end + ###################################################################### # # # Multi-PANEL INDUCED VELOCITIES # @@ -432,6 +848,14 @@ end @testset "Multi-Panel, Multi-Target Induced Velocity Matrices" begin @testset "Multiple Vortex Panel Induced Velocities on Multiple Points" begin + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + + # define control points controlpoints = [0.0 1.0; 1.0 1.0]' @@ -447,26 +871,52 @@ end # use unit strengths strengths = [1.0 1.0; 1.0 1.0] + # - GaussKronrod - # # [cp, n, x/r] VEL = dt.induced_velocities_from_vortex_panels_on_points( - controlpoints, nodes, nodemap, influence_lengths, strengths + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=gk_integration_options, singular=gk_integration_options), ) # [vz1 vr1; vz2 vr2] vn12cp1 = dt.self_vortex_panel_integration( - nodes[:, 1], nodes[:, 2], influence_lengths[1], controlpoints[:, 1], zeros(20) + gk_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + gk_cache, ) vn12cp2 = dt.nominal_vortex_panel_integration( - nodes[:, 1], nodes[:, 2], influence_lengths[1], controlpoints[:, 2], zeros(20) + gk_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + gk_cache, ) vn23cp1 = dt.nominal_vortex_panel_integration( - nodes[:, 2], nodes[:, 3], influence_lengths[1], controlpoints[:, 1], zeros(20) + gk_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + gk_cache, ) vn23cp2 = dt.nominal_vortex_panel_integration( - nodes[:, 2], nodes[:, 3], influence_lengths[1], controlpoints[:, 2], zeros(20) + gk_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + gk_cache, ) @test all(VEL[1, 1, :] .== vn12cp1[1, :]) @@ -476,10 +926,131 @@ end @test all(VEL[2, 1, :] .== vn12cp2[1, :]) @test all(VEL[2, 2, :] .== vn12cp2[2, :] .+ vn23cp2[1, :]) @test all(VEL[2, 3, :] .== vn23cp2[2, :]) + + # - GaussLegendre - # + # [cp, n, x/r] + VEL = dt.induced_velocities_from_vortex_panels_on_points( + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=gl_integration_options, singular=gl_integration_options), + ) + + # [vz1 vr1; vz2 vr2] + vn12cp1 = dt.self_vortex_panel_integration( + gl_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + gl_cache, + ) + + vn12cp2 = dt.nominal_vortex_panel_integration( + gl_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + gl_cache, + ) + + vn23cp1 = dt.nominal_vortex_panel_integration( + gl_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + gl_cache, + ) + + vn23cp2 = dt.nominal_vortex_panel_integration( + gl_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + gl_cache, + ) + + @test all(VEL[1, 1, :] .== vn12cp1[1, :]) + #note that node connected to 2 panels has influence contributions from both panels + @test all(VEL[1, 2, :] .== vn12cp1[2, :] .+ vn23cp1[1, :]) + @test all(VEL[1, 3, :] .== vn23cp1[2, :]) + @test all(VEL[2, 1, :] .== vn12cp2[1, :]) + @test all(VEL[2, 2, :] .== vn12cp2[2, :] .+ vn23cp2[1, :]) + @test all(VEL[2, 3, :] .== vn23cp2[2, :]) + + # - Romberg - # + # [cp, n, x/r] + VEL = dt.induced_velocities_from_vortex_panels_on_points( + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=r_integration_options, singular=r_integration_options), + ) + + # [vz1 vr1; vz2 vr2] + vn12cp1 = dt.self_vortex_panel_integration( + r_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + r_cache, + ) + + vn12cp2 = dt.nominal_vortex_panel_integration( + r_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + r_cache, + ) + + vn23cp1 = dt.nominal_vortex_panel_integration( + r_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + r_cache, + ) + + vn23cp2 = dt.nominal_vortex_panel_integration( + r_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + r_cache, + ) + + @test all(VEL[1, 1, :] .== vn12cp1[1, :]) + #note that node connected to 2 panels has influence contributions from both panels + @test all(VEL[1, 2, :] .== vn12cp1[2, :] .+ vn23cp1[1, :]) + @test all(VEL[1, 3, :] .== vn23cp1[2, :]) + #TODO: why is the answer different when directly vs indirectly calling integration? + @test all(isapprox.(VEL[2, 1, :] , vn12cp2[1, :], atol=1e-4)) + @test all(isapprox.(VEL[2, 2, :] , vn12cp2[2, :] .+ vn23cp2[1, :], atol=1e-4)) + @test all(VEL[2, 3, :] .== vn23cp2[2, :]) end @testset "Multiple Source Panel Induced Velocities on Multiple Points" begin + gk_integration_options = dt.GaussKronrod() + gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) + gl_integration_options = dt.GaussLegendre() + gl_cache = dt.allocate_integration_containers(gl_integration_options, 1.0) + r_integration_options = dt.Romberg() + r_cache = dt.allocate_integration_containers(r_integration_options, 1.0) + + # define control points controlpoints = [0.0 1.0; 1.0 1.0]' @@ -495,26 +1066,107 @@ end # use unit strengths strengths = [1.0 1.0; 1.0 1.0] + # GaussKronrod + # [cp, n, x/r] + VEL = dt.induced_velocities_from_source_panels_on_points( + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=gk_integration_options, singular=gk_integration_options), + ) + + # [vz1 vr1; vz2 vr2] + vn12cp1 = dt.self_source_panel_integration( + gk_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + gk_cache, + ) + + vn12cp2 = dt.nominal_source_panel_integration( + gk_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + gk_cache, + ) + + vn23cp1 = dt.nominal_source_panel_integration( + gk_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + gk_cache, + ) + + vn23cp2 = dt.nominal_source_panel_integration( + gk_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + gk_cache, + ) + + @test all(VEL[1, 1, :] .== vn12cp1[1, :]) + @test all(VEL[1, 2, :] .== vn12cp1[2, :] .+ vn23cp1[1, :]) + @test all(VEL[1, 3, :] .== vn23cp1[2, :]) + @test all(VEL[2, 1, :] .== vn12cp2[1, :]) + @test all(VEL[2, 2, :] .== vn12cp2[2, :] .+ vn23cp2[1, :]) + @test all(VEL[2, 3, :] .== vn23cp2[2, :]) + + # GaussLegendre # [cp, n, x/r] VEL = dt.induced_velocities_from_source_panels_on_points( - controlpoints, nodes, nodemap, influence_lengths, strengths + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=gl_integration_options, singular=gl_integration_options), ) # [vz1 vr1; vz2 vr2] vn12cp1 = dt.self_source_panel_integration( - nodes[:, 1], nodes[:, 2], influence_lengths[1], controlpoints[:, 1], zeros(20) + gl_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + gl_cache, ) vn12cp2 = dt.nominal_source_panel_integration( - nodes[:, 1], nodes[:, 2], influence_lengths[1], controlpoints[:, 2], zeros(20) + gl_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + gl_cache, ) vn23cp1 = dt.nominal_source_panel_integration( - nodes[:, 2], nodes[:, 3], influence_lengths[1], controlpoints[:, 1], zeros(20) + gl_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + gl_cache, ) vn23cp2 = dt.nominal_source_panel_integration( - nodes[:, 2], nodes[:, 3], influence_lengths[1], controlpoints[:, 2], zeros(20) + gl_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + gl_cache, ) @test all(VEL[1, 1, :] .== vn12cp1[1, :]) @@ -523,5 +1175,61 @@ end @test all(VEL[2, 1, :] .== vn12cp2[1, :]) @test all(VEL[2, 2, :] .== vn12cp2[2, :] .+ vn23cp2[1, :]) @test all(VEL[2, 3, :] .== vn23cp2[2, :]) + + # Romberg + # [cp, n, x/r] + VEL = dt.induced_velocities_from_source_panels_on_points( + controlpoints, + nodes, + nodemap, + influence_lengths, + strengths, + (; nominal=r_integration_options, singular=r_integration_options), + ) + + # [vz1 vr1; vz2 vr2] + vn12cp1 = dt.self_source_panel_integration( + r_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 1], + r_cache, + ) + + vn12cp2 = dt.nominal_source_panel_integration( + r_integration_options, + nodes[:, 1], + nodes[:, 2], + influence_lengths[1], + controlpoints[:, 2], + r_cache, + ) + + vn23cp1 = dt.nominal_source_panel_integration( + r_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 1], + r_cache, + ) + + vn23cp2 = dt.nominal_source_panel_integration( + r_integration_options, + nodes[:, 2], + nodes[:, 3], + influence_lengths[1], + controlpoints[:, 2], + r_cache, + ) + + @test all(VEL[1, 1, :] .== vn12cp1[1, :]) + @test all(VEL[1, 2, :] .== vn12cp1[2, :] .+ vn23cp1[1, :]) + @test all(VEL[1, 3, :] .== vn23cp1[2, :]) + # TODO: why are these specific two different (but not totally wrong)?? + @test all(isapprox.(VEL[2, 1, :] , vn12cp2[1, :], atol=1e-4)) + @test all(isapprox.(VEL[2, 2, :] , vn12cp2[2, :] .+ vn23cp2[1, :],atol=1e-3)) + @test all(VEL[2, 3, :] .== vn23cp2[2, :]) end end diff --git a/test/linear_system_assembly.jl b/test/linear_system_assembly.jl index a175565d..42b52a5d 100644 --- a/test/linear_system_assembly.jl +++ b/test/linear_system_assembly.jl @@ -1,5 +1,6 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") @testset "Linear System Assembly" begin + integration_options = (; nominal=dt.GaussLegendre(), singular=dt.GaussLegendre()) # define coordinates x1 = [1.0; 0.5; 0.0; 0.5; 1.0] @@ -25,6 +26,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.node, panels.nodemap, panels.influence_length, + integration_options, ) dt.add_te_gap_aic!( @@ -36,6 +38,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.tendotn, panels.tencrossn, panels.teadjnodeidxs, + integration_options ) AICpcp = dt.vortex_aic_boundary_on_field( @@ -44,6 +47,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.node, panels.nodemap, panels.influence_length, + integration_options ) dt.add_te_gap_aic!( @@ -55,6 +59,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.tendotn, panels.tencrossn, panels.teadjnodeidxs, + integration_options ) Vinf = 1.0 #magnitude doesn't matter yet. @@ -129,6 +134,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.node, panels.nodemap, panels.influence_length, + integration_options ) AICpcp = dt.vortex_aic_boundary_on_field( @@ -137,6 +143,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") panels.node, panels.nodemap, panels.influence_length, + integration_options ) Vinf = 1.0 #magnitude doesn't matter yet. diff --git a/test/runtests.jl b/test/runtests.jl index 5e1aac8f..2f724694 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,16 +23,13 @@ include("test_utils.jl") println("Running Tests...") -#---------------------------------# -# New Tests # -#---------------------------------# - # - Caching Tests - # include("iotests.jl") # - pre-process related tests - # include("afcorrections.jl") include("panel_generation_tests.jl") +include("induced_velocities.jl") include("influence_coefficients.jl") include("linear_system_assembly.jl") include("pre_processing_tests.jl")