update code and tests to allow for various implemented quadrature opt…

…ions
byuflowlab · Apr 4, 2024 · 9a81bfa · 9a81bfa
1 parent c97f772
commit 9a81bfa
Show file tree

Hide file tree

Showing 11 changed files with 795 additions and 75 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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