diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index 29350e7b..8c68c0ff 100644 --- a/src/DuctAPE.jl +++ b/src/DuctAPE.jl @@ -80,11 +80,13 @@ include("utilities/airfoils/airfoil_utilities.jl") include("utilities/airfoils/naca_65series.jl") ##### ----- Analysis ----- ##### +include("analysis/setup.jl") include("analysis/analyses.jl") ##### ----- PREPROCESS ----- ##### # Pre-solve initializations include("preprocess/preprocess.jl") +include("preprocess/initialize_states.jl") # Geometry Functions include("preprocess/geometry/body_geometry.jl") @@ -106,6 +108,7 @@ include("preprocess/velocities/gausskronrod_integrals.jl") ##### ----- PROCESS ----- ##### # Solve and Residual Functions +include("process/process.jl") include("process/solve.jl") include("process/residuals/CSORresidual.jl") include("process/residuals/systemresidual.jl") diff --git a/src/analysis/analyses.jl b/src/analysis/analyses.jl index 6e6ec988..f5401f41 100644 --- a/src/analysis/analyses.jl +++ b/src/analysis/analyses.jl @@ -2,104 +2,20 @@ """ function analyze( propulsor::Propulsor, - options=set_options(); - # precomp_container_caching=nothing, + options::Options=set_options(); + prepost_container_caching=nothing, solve_parameter_caching=nothing, solve_container_caching=nothing, - # post_caching=nothing + return_inputs=false, ) - # - Get type to dispatch caches - # - TF = promote_type( - eltype(propulsor.duct_coordinates), - eltype(propulsor.centerbody_coordinates), - eltype(propulsor.operating_point.Vinf), - eltype(propulsor.operating_point.Omega), - eltype(propulsor.operating_point.rhoinf), - eltype(propulsor.operating_point.muinf), - eltype(propulsor.operating_point.asound), - eltype(propulsor.rotorstator_parameters.B), - eltype(propulsor.rotorstator_parameters.Rhub), - eltype(propulsor.rotorstator_parameters.Rtip), - eltype(propulsor.rotorstator_parameters.rotorzloc), - eltype(propulsor.rotorstator_parameters.chords), - eltype(propulsor.rotorstator_parameters.twists), - ) - - # if isnothing(precomp_container_caching) - # precomp_container_caching = allocate_precomp_container_cache( - # propulsor.paneling_constants - # ) - # end - - # if isnothing(post_caching) - # post_caching = allocate_post_cache( - # propulsor.paneling_constants - # ) - # end - - # - Pull out the Caches - # - # (; precomp_container_cache, precomp_container_cache_dims) = precomp_container_caching - # (; post_cache, post_cache_dims) = post_caching - - # # - Reshape precomp_container_cache - # - # precomp_container_cache_vec = @views PreallocationTools.get_tmp(precomp_container_cache, TF(1.0)) - # precomp_containers = withdraw_precomp_container_cache( - # precomp_container_cache, precomp_container_cache_dims - # ) - - # - Set up Solver Sensitivity Paramter Cache - # - - # Allocate Cache - if isnothing(solve_parameter_caching) - solve_parameter_caching = allocate_solve_parameter_cache( - options.solver_options, propulsor.paneling_constants - ) - end - - # separate out caching items - (; solve_parameter_cache, solve_parameter_cache_dims) = solve_parameter_caching - - # get correct cache type - solve_parameter_cache_vector = @views PreallocationTools.get_tmp( - solve_parameter_cache, TF(1.0) - ) - # reset cache - solve_parameter_cache_vector .= 0 - - # reshape cache - solve_parameter_tuple = withdraw_solve_parameter_cache( - options.solver_options, solve_parameter_cache_vector, solve_parameter_cache_dims - ) - - # copy over operating point - for f in fieldnames(typeof(propulsor.operating_point)) - solve_parameter_tuple.operating_point[f] .= getfield(propulsor.operating_point, f) - end - - # - Do precomputations - # - if options.verbose - println("Pre-computing Parameters") - end - - # 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!( - 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.grid_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, + # - Set Up - # + problem_dimensions, prepost_containers, solve_parameter_cache_vector, solve_parameter_cache_dims, ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps = setup_analysis( + propulsor, + options; + prepost_container_caching=prepost_container_caching, + solve_parameter_caching=solve_parameter_caching, + solve_container_caching=solve_container_caching, ) # - Check that the precomputation went well - # @@ -115,6 +31,7 @@ function analyze( end end #TODO: write a function that returns the same as outs below, but all zeros + #TODO: probably just call the post-process function directly and return a reset_container! of the output return [],#zero_outputs(), (; solve_parameter_tuple..., ivb, airfoils, idmaps, panels, problem_dimensions), false @@ -123,17 +40,17 @@ function analyze( # - Continue with Analysis - # return analyze( propulsor, + prepost_containers, solve_parameter_cache_vector, solve_parameter_cache_dims, airfoils, ivb, A_bb_LU, - panels, idmaps, problem_dimensions, options; solve_container_caching=solve_container_caching, - # post_caching=nothing + return_inputs=return_inputs, ) end @@ -141,23 +58,19 @@ end """ function analyze( propulsor::Propulsor, + prepost_containers, solve_parameter_cache_vector, solve_parameter_cache_dims, airfoils, ivb, A_bb_LU, - panels, idmaps, problem_dimensions, - options=set_options(); + options::Options=set_options(); return_inputs=false, - # precomp_container_caching=nothing, solve_container_caching=nothing, - # post_caching=nothing ) - # - Finish Pre-Processing - # - # Set up Solve Container Cache if isnothing(solve_container_caching) solve_container_caching = allocate_solve_container_cache( @@ -180,16 +93,15 @@ function analyze( # - Post-Process - # outs = post_process( options.solver_options, - solve_container_caching, converged_states, + prepost_containers, + solve_container_caching, solve_parameter_cache_vector, solve_parameter_cache_dims, propulsor.operating_point, propulsor.reference_parameters, - ivb, A_bb_LU, airfoils, - panels, idmaps, problem_dimensions; write_outputs=options.write_outputs, @@ -206,6 +118,8 @@ function analyze( return outs, (; + prepost_containers.panels, + prepost_containers.ivb, solve_parameter_tuple..., blade_elements=(; blade_elements..., airfoils...), linsys=(; linsys..., A_bb_LU), @@ -218,135 +132,195 @@ end """ """ -function process( - solver_options::TS, - solve_parameter_cache_vector, - solve_parameter_cache_dims, - airfoils, - A_bb_LU, - solve_container_caching, - idmaps, - options, -) where {TS<:Union{ExternalSolverOptions,MultiSolverOptions}} +function analyze( + multipoint::AbstractVector{TO}, + propulsor::Propulsor, + options::Options=set_options(); + prepost_container_caching=nothing, + solve_parameter_caching=nothing, + solve_container_caching=nothing, + return_inputs=false, +) where {TO<:OperatingPoint} - # - Initialize Aero - # - if options.verbose - println("\nInitializing Velocities") - end - # view the initial conditions out of the inputs cache - solve_parameter_tuple = withdraw_solve_parameter_cache( - solver_options, solve_parameter_cache_vector, solve_parameter_cache_dims + # - Set Up - # + problem_dimensions, prepost_containers, solve_parameter_cache_vector, solve_parameter_cache_dims, ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps = setup_analysis( + propulsor, + options; + prepost_container_caching=prepost_container_caching, + solve_parameter_caching=solve_parameter_caching, + solve_container_caching=solve_container_caching, ) - (; vz_rotor, vtheta_rotor, Cm_wake, operating_point, linsys, ivr, ivw) = - solve_parameter_tuple - # initialize velocities - # TODO; add some sort of unit test for this function - initialize_velocities!( - vz_rotor, - vtheta_rotor, - Cm_wake, - solve_parameter_tuple.operating_point, - (; solve_parameter_tuple.blade_elements..., airfoils...), - (; solve_parameter_tuple.linsys..., A_bb_LU), - solve_parameter_tuple.ivr, - solve_parameter_tuple.ivw, - idmaps.body_totnodes, - idmaps.wake_panel_sheet_be_map, + # - Check that the precomputation went well - # + #= + 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.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.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 + #TODO: write a function that returns the same as outs below, but all zeros + #TODO: probably just call the post-process function directly and return a reset_container! of the output + return [],#zero_outputs(), + (; solve_parameter_tuple..., ivb, airfoils, idmaps, panels, problem_dimensions), + false + end + + # - initialize converged state multi-point vector - # + converged_states = zeros( + promote_propulosor_type(propulsor), + solve_parameter_cache_dims.state_dims[end].index[end], + length(multipoint), ) - # - combine cache and constants - # - const_cache = (; - # - Constants - # - options.verbose, - options.silence_warnings, - #nlsolve options - solver_options, - # Constant Parameters + return analyze( + multipoint, + propulsor, + prepost_containers, + solve_parameter_cache_vector, + solve_parameter_cache_dims, airfoils, + ivb, A_bb_LU, idmaps, - solve_parameter_cache_dims, - # Cache(s) - solve_container_caching..., - solve_parameter_tuple..., # contains SIAMFANLEOptions containers needed for solver definition - ) - - # - Solve with ImplicitAD - # - if options.verbose - println("\nSolving Nonlinear System") - end - return ImplicitAD.implicit( - solve, system_residual!, solve_parameter_cache_vector, const_cache + problem_dimensions, + options, + converged_states; + solve_container_caching=solve_container_caching, + return_inputs=return_inputs, ) end -""" -""" -function process( - solver_options::CSORSolverOptions, +function analyze( + multipoint::AbstractVector{TO}, + propulsor::Propulsor, + prepost_containers, solve_parameter_cache_vector, solve_parameter_cache_dims, airfoils, + ivb, A_bb_LU, - solve_container_caching, idmaps, - options, -) - - # - Initialize Aero - # + problem_dimensions, + options::Options, + converged_states; + solve_container_caching=nothing, + return_inputs=false, +) where {TO<:OperatingPoint} if options.verbose - println("\nInitializing Velocities") + println("\nRunning Multipoint Analysis") + end + + # Set up Solve Container Cache + if isnothing(solve_container_caching) + solve_container_caching = allocate_solve_container_cache( + options.solver_options, propulsor.paneling_constants + ) end - # view the initial conditions out of the inputs cache + solve_parameter_tuple = withdraw_solve_parameter_cache( - solver_options, solve_parameter_cache_vector, solve_parameter_cache_dims + options.solver_options, solve_parameter_cache_vector, solve_parameter_cache_dims ) - (; Gamr, sigr, gamw, operating_point, blade_elements, linsys, ivr, ivw, wakeK) = - solve_parameter_tuple - # - Initialize States - # - initialize_strengths!( - Gamr, - sigr, - gamw, - operating_point, - (; blade_elements..., airfoils...), - (; linsys..., A_bb_LU), - ivr, - ivw, - wakeK, - idmaps.body_totnodes, - idmaps.wake_nodemap, - idmaps.wake_endnodeidxs, - idmaps.wake_panel_sheet_be_map, - idmaps.wake_node_sheet_be_map, - idmaps.wake_node_ids_along_casing_wake_interface, - idmaps.wake_node_ids_along_centerbody_wake_interface, - ) + # reset multipoint index + options.multipoint_index[] = 0 - # - combine cache and constants - # - const_cache = (; - # - Constants - # - options.silence_warnings, - options.verbose, - #CSOR solve options - solver_options, - # Constant Parameters - airfoils, - A_bb_LU, - idmaps, - # Cache(s) - solve_parameter_cache_vector, - solve_parameter_cache_dims, - solve_container_caching..., - ) + # - Loop through Operating Points - # + for (iop, op) in enumerate(multipoint) - # - Solve with ImplicitAD - # - if options.verbose - println("\nSolving Nonlinear System using CSOR Method") + # incremenet operating point to keep track of convergence in multipoint solves + options.multipoint_index[] += 1 + + if options.verbose + println("\n Operating Point:") + for fn in fieldnames(typeof(op)) + println(@sprintf " %6s = %5.3e" fn getfield(op, fn)[]) + # println(" ", fn, " = ", getfield(op,fn)) + end + end + + # - copy over operating point - # + for f in fieldnames(typeof(op)) + solve_parameter_tuple.operating_point[f] .= getfield(op, f) + end + + # - Set up Body Linear System RHS - # + vinfvec = [op.Vinf[]; 0.0] + prepost_containers.vdnb .= [ + dot(vinfvec, nhat) for + nhat in eachcol(prepost_containers.panels.body_vortex_panels.normal) + ] + prepost_containers.vdnpcp .= [ + dot(vinfvec, nhat) for + nhat in eachcol(prepost_containers.panels.body_vortex_panels.itnormal) + ] + assemble_rhs_matrix!( + solve_parameter_tuple.linsys.b_bf, + prepost_containers.vdnb, + prepost_containers.vdnpcp, + prepost_containers.panels.body_vortex_panels.npanel, + prepost_containers.panels.body_vortex_panels.nnode, + prepost_containers.panels.body_vortex_panels.totpanel, + prepost_containers.panels.body_vortex_panels.totnode, + prepost_containers.panels.body_vortex_panels.prescribednodeidxs, + ) + + # - Process - # + converged_states[:, iop] .= process( + options.solver_options, + solve_parameter_cache_vector, + solve_parameter_cache_dims, + airfoils, + A_bb_LU, + solve_container_caching, + idmaps, + options, + ) + end + + # - Post-Process - # + outs = [ + post_process( + options.solver_options, + cs, + prepost_containers, + solve_container_caching, + solve_parameter_cache_vector, + solve_parameter_cache_dims, + propulsor.operating_point, + propulsor.reference_parameters, + A_bb_LU, + airfoils, + idmaps, + problem_dimensions; + write_outputs=options.write_outputs[iop], + outfile=options.outfile[iop], + checkoutfileexists=options.checkoutfileexists, + output_tuple_name=options.output_tuple_name[iop], + verbose=options.verbose, + ) for (iop, cs) in enumerate(eachcol(converged_states)) + ] + + if return_inputs + solve_parameter_tuple = withdraw_solve_parameter_cache( + solver_options, solve_parameter_cache_vector, solve_parameter_cache_dims + ) + + return outs, + (; + prepost_containers.panels, + prepost_containers.ivb, + solve_parameter_tuple..., + blade_elements=(; blade_elements..., airfoils...), + linsys=(; linsys..., A_bb_LU), + ), + options.solver_options.converged + else + return outs, options.solver_options.converged end - return ImplicitAD.implicit( - solve, CSOR_residual!, solve_parameter_cache_vector, const_cache - ) end + diff --git a/src/analysis/setup.jl b/src/analysis/setup.jl new file mode 100644 index 00000000..5f971f8b --- /dev/null +++ b/src/analysis/setup.jl @@ -0,0 +1,122 @@ +""" +""" +function setup_analysis( + propulsor::Propulsor, + options=set_options(); + prepost_container_caching=nothing, + solve_parameter_caching=nothing, + solve_container_caching=nothing, +) + # - Get type to dispatch caches - # + TF = promote_type( + eltype(propulsor.duct_coordinates), + eltype(propulsor.centerbody_coordinates), + eltype(propulsor.operating_point.Vinf), + eltype(propulsor.operating_point.Omega), + eltype(propulsor.operating_point.rhoinf), + eltype(propulsor.operating_point.muinf), + eltype(propulsor.operating_point.asound), + eltype(propulsor.rotorstator_parameters.B), + eltype(propulsor.rotorstator_parameters.Rhub), + eltype(propulsor.rotorstator_parameters.Rtip), + eltype(propulsor.rotorstator_parameters.rotorzloc), + eltype(propulsor.rotorstator_parameters.chords), + eltype(propulsor.rotorstator_parameters.twists), + ) + + # - Get Problem Dimensions - # + problem_dimensions = get_problem_dimensions(propulsor.paneling_constants) + + ##### ----- SET UP CACHES AS NEEDED ----- ##### + + # - Set up Pre- and Post-process Cache - # + # Allocate Cache + if isnothing(prepost_container_caching) + prepost_container_caching = allocate_prepost_container_cache( + propulsor.paneling_constants + ) + end + + # unpack the caching + (; prepost_container_cache, prepost_container_cache_dims) = prepost_container_caching + + # Get correct cached types + prepost_container_cache_vec = @views PreallocationTools.get_tmp( + prepost_container_cache, TF(1.0) + ) + + # reset cache + prepost_container_cache_vec .= 0 + + # Reshape Cache + prepost_containers = withdraw_prepost_container_cache( + prepost_container_cache_vec, prepost_container_cache_dims + ) + + # - Set up Solver Sensitivity Paramter Cache - # + + # Allocate Cache + if isnothing(solve_parameter_caching) + solve_parameter_caching = allocate_solve_parameter_cache( + options.solver_options, propulsor.paneling_constants + ) + end + + # unpack caching + (; solve_parameter_cache, solve_parameter_cache_dims) = solve_parameter_caching + + # get correct cache type + solve_parameter_cache_vector = @views PreallocationTools.get_tmp( + solve_parameter_cache, TF(1.0) + ) + # reset cache + solve_parameter_cache_vector .= 0 + + # reshape cache + solve_parameter_tuple = withdraw_solve_parameter_cache( + options.solver_options, solve_parameter_cache_vector, solve_parameter_cache_dims + ) + + # copy over operating point + for f in fieldnames(typeof(propulsor.operating_point)) + solve_parameter_tuple.operating_point[f] .= getfield(propulsor.operating_point, f) + end + + # - Do preprocessutations - # + if options.verbose + println("Pre-computing Parameters") + end + + ##### ----- PERFORM PREPROCESSING COMPUTATIONS ----- ##### + + # - Preprocess - # + ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps, _ = precompute_parameters!( + solve_parameter_tuple.ivr, + solve_parameter_tuple.ivw, + solve_parameter_tuple.blade_elements, + solve_parameter_tuple.linsys, + solve_parameter_tuple.wakeK, + propulsor, + prepost_containers, + problem_dimensions; + grid_solver_options=options.grid_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, + ) + + return problem_dimensions, + prepost_containers, + solve_parameter_cache_vector, + solve_parameter_cache_dims, + ivb, + A_bb_LU, + lu_decomp_flag, + airfoils, + idmaps +end diff --git a/src/postprocess/postprocess.jl b/src/postprocess/postprocess.jl index 987bb08f..77e8f4bf 100644 --- a/src/postprocess/postprocess.jl +++ b/src/postprocess/postprocess.jl @@ -1,134 +1,23 @@ -""" -""" -function run_residual!( - solver_options::TS, - state_variables, - state_dims, - solve_container_cache, - solve_container_cache_dims, - operating_point, - ivr, - ivw, - linsys, - blade_elements, - wakeK, - idmaps, -) where {TS<:ExternalSolverOptions} - - #= - NOTE: we want to get all the intermediate values available to user if desired. - The solve_containers cache will contain all the intermediate values after running the estimate states function. - =# - # - Separate out the state variables - # - vz_rotor, vtheta_rotor, Cm_wake = extract_state_variables( - solver_options, state_variables, state_dims - ) - - # - Extract and Reset Cache - # - # get cache vector of correct types - solve_container_cache_vec = @views PreallocationTools.get_tmp( - solve_container_cache, state_variables - ) - solve_containers = withdraw_solve_container_cache( - solver_options, solve_container_cache_vec, solve_container_cache_dims - ) - reset_containers!(solve_containers) #note: also zeros out state estimates - - # - Estimate New States - # - # currently has 280 allocations - estimate_states!( - solve_containers, - vz_rotor, - vtheta_rotor, - Cm_wake, - operating_point, - ivr, - ivw, - linsys, - blade_elements, - wakeK, - idmaps, - ) - - return (; vz_rotor, vtheta_rotor, Cm_wake, solve_containers...) -end - -""" -""" -function run_residual!( - solver_options::CSORSolverOptions, - state_variables, - state_dims, - solve_container_cache, - solve_container_cache_dims, - operating_point, - ivr, - ivw, - linsys, - blade_elements, - wakeK, - idmaps, -) - - #= - NOTE: we want to get all the intermediate values available to user if desired. - The solve_containers cache will contain all the intermediate values after running the insides of the residual function - =# - # - Separate out the state variables - # - Gamr, sigr, gamw = extract_state_variables(solver_options, state_variables, state_dims) - - # - Extract and Reset Cache - # - # get cache vector of correct types - solve_container_cache_vec = @views PreallocationTools.get_tmp( - solve_container_cache, state_variables - ) - solve_container_cache_vec .= 0 - solve_containers = withdraw_solve_container_cache( - solver_options, solve_container_cache_vec, solve_container_cache_dims - ) - - # - Run Residual - # - compute_CSOR_residual!( - zeros(2), - solver_options, - solve_containers, - Gamr, - sigr, - gamw, - operating_point, - ivr, - ivw, - linsys, - blade_elements, - wakeK, - idmaps; - verbose=false, - ) - - return (; Gamr, sigr, gamw, solve_containers...) -end - """ """ function post_process( solver_options, + converged_states, + prepost_containers, solve_container_caching, - state_variables, solve_parameter_cache_vector, solve_parameter_cache_dims, operating_point, reference_parameters, - ivb, A_bb_LU, airfoils, - panels, idmaps, problem_dimensions; - write_outputs=false, - outfile="outputs.jl", - checkoutfileexists=false, - output_tuple_name="outs", - verbose=false, + write_outputs=options.write_outputs, + outfile=options.outfile, + checkoutfileexists=options.checkoutfileexists, + output_tuple_name=options.output_tuple_name, + verbose=options.verbose, ) ### --- SETUP --- ### @@ -139,6 +28,54 @@ function post_process( # - Extract Reference Parameters - # (; Vref, Rref) = reference_parameters + # - Extract PrePost Cache - # + reset_containers!(prepost_containers; exception_keys=(:panels,:ivb)) + (; + # stuff from pre-process + panels, + ivb, + # rotor stuff + rotor_inviscid_thrust, + rotor_inviscid_thrust_dist, + rotor_viscous_thrust, + rotor_viscous_thrust_dist, + rotor_thrust, + rotor_inviscid_torque, + rotor_inviscid_torque_dist, + rotor_viscous_torque, + rotor_viscous_torque_dist, + rotor_torque, + rotor_inviscid_power, + rotor_inviscid_power_dist, + rotor_viscous_power, + rotor_viscous_power_dist, + rotor_power, + rotor_CT, + rotor_CQ, + rotor_CP, + rotor_efficiency, + induced_efficiency, + blade_normal_force_per_unit_span, + blade_tangential_force_per_unit_span, + blade_loading_intermediate_containers, + # body stuff + zpts, + vtan_tuple, + cp_tuple, + body_thrust, + body_force_coefficient, + # cp_tuple, + # totals stuff + total_thrust, + total_torque, + total_power, + total_efficiency, + ideal_efficiency, + total_CT, + total_CQ, + total_CP, + ) = prepost_containers + # - Extract Panels - # (; body_vortex_panels, rotor_source_panels, wake_vortex_panels) = panels @@ -159,7 +96,6 @@ function post_process( (; solve_container_cache, solve_container_cache_dims) = solve_container_caching # - Run Residual to get intermediate values - # - # for combination solvers, get solver type that was last run if typeof(solver_options) <: CompositeSolverOptions idx = findlast(x -> x, (p -> p.converged[1]).(solver_options.solvers)) @@ -173,7 +109,7 @@ function post_process( res_vals = run_residual!( sopt, - state_variables, + converged_states, solve_parameter_cache_dims.state_dims, solve_container_cache, solve_container_cache_dims, @@ -207,33 +143,53 @@ function post_process( ### ----- ROTOR OUTPUTS ----- ### # rename for convenience (; nbe, nrotor) = problem_dimensions - rotor_panel_lengths = reshape(rotor_panel_lengths, (nbe, nrotor)) rotor_panel_centers = blade_elements.rotor_panel_centers B = blade_elements.B chords = blade_elements.chords # - Rotor Thrust - # # inviscid thrust - rotor_inviscid_thrust, rotor_inviscid_thrust_dist = inviscid_rotor_thrust( - Ctheta_rotor, Gamma_tilde, rotor_panel_lengths, rhoinf[1] + inviscid_rotor_thrust!( + rotor_inviscid_thrust, + rotor_inviscid_thrust_dist, + Ctheta_rotor, + Gamma_tilde, + rotor_panel_lengths, + rhoinf[1], ) # viscous thrust - rotor_viscous_thrust, rotor_viscous_thrust_dist = viscous_rotor_thrust( - Cz_rotor, Cmag_rotor, B, chords, rotor_panel_lengths, cd, rhoinf[1] + viscous_rotor_thrust!( + rotor_viscous_thrust, + rotor_viscous_thrust_dist, + Cz_rotor, + Cmag_rotor, + B, + chords, + rotor_panel_lengths, + cd, + rhoinf[1], ) # total thrust - rotor_thrust = rotor_inviscid_thrust .+ rotor_viscous_thrust + rotor_thrust .= rotor_inviscid_thrust .+ rotor_viscous_thrust # - Rotor Torque - # # inviscid torque - rotor_inviscid_torque, rotor_inviscid_torque_dist = inviscid_rotor_torque( - Cz_rotor, rotor_panel_centers, rotor_panel_lengths, Gamma_tilde, rhoinf[1] + inviscid_rotor_torque!( + rotor_inviscid_torque, + rotor_inviscid_torque_dist, + Cz_rotor, + rotor_panel_centers, + rotor_panel_lengths, + Gamma_tilde, + rhoinf[1], ) # viscous torque - rotor_viscous_torque, rotor_viscous_torque_dist = viscous_rotor_torque( + viscous_rotor_torque!( + rotor_viscous_torque, + rotor_viscous_torque_dist, Ctheta_rotor, Cmag_rotor, B, @@ -245,64 +201,81 @@ function post_process( ) # total torque - rotor_torque = rotor_inviscid_torque .+ rotor_viscous_torque + rotor_torque .= rotor_inviscid_torque .+ rotor_viscous_torque # - Rotor Power - # # inviscid power - rotor_inviscid_power = inviscid_rotor_power(rotor_inviscid_torque', Omega) - - rotor_inviscid_power_dist = similar(rotor_inviscid_torque_dist) .= 0.0 - for ir in 1:nrotor - rotor_inviscid_power_dist[:, ir] = inviscid_rotor_power( - @view(rotor_inviscid_torque_dist[:, ir]), Omega[ir] - ) - end + rotor_power!( + rotor_inviscid_power, + rotor_inviscid_power_dist, + rotor_inviscid_torque, + rotor_inviscid_torque_dist, + Omega, + ) # viscous power - rotor_viscous_power = viscous_rotor_power(rotor_viscous_torque', Omega) - rotor_viscous_power_dist = similar(rotor_viscous_torque_dist) .= 0.0 - for ir in 1:nrotor - rotor_inviscid_power_dist[:, ir] = viscous_rotor_power( - @view(rotor_inviscid_torque_dist[:, ir]), Omega[ir] - ) - end + rotor_power!( + rotor_viscous_power, + rotor_viscous_power_dist, + rotor_viscous_torque, + rotor_viscous_torque_dist, + Omega, + ) # total power - rotor_power = rotor_inviscid_power .+ rotor_viscous_power + rotor_power .= rotor_inviscid_power .+ rotor_viscous_power # - Rotor Performance Coefficients - # - rotor_CT, rotor_CQ, rotor_CP = tqpcoeff( - rotor_thrust, rotor_torque, rotor_power, rhoinf[1], Omega, Rref[1] + tqpcoeff!( + rotor_CT, + rotor_CQ, + rotor_CP, + rotor_thrust, + rotor_torque, + rotor_power, + rhoinf[1], + Omega, + Rref[1], ) # - Rotor Efficiency - # - rotor_efficiency = get_total_efficiency(rotor_thrust, rotor_power, Vinf[1]) + get_total_efficiency!(rotor_efficiency, rotor_thrust, rotor_power, Vinf[1]) # - Blade Loading - # - blade_normal_force_per_unit_span, blade_tangential_force_per_unit_span = get_blade_loads( - Cmag_rotor, beta1, cl, cd, chords, rhoinf[1] + get_blade_loads!( + blade_normal_force_per_unit_span, + blade_tangential_force_per_unit_span, + Cmag_rotor, + beta1, + cl, + cd, + chords, + rhoinf[1], + blade_loading_intermediate_containers, ) ### --- BODY OUTPUTS --- ### # - Surface Velocity on Bodies - # # TODO: update any tests for this function - vtan_tuple = get_body_tangential_velocities( + get_body_tangential_velocities!( + vtan_tuple, gamb, gamw, sigr, ivb, - Vinf[1], - body_vortex_panels.totnode, - body_vortex_panels.totpanel, - body_vortex_panels.nnode, - body_vortex_panels.npanel, + Vinf[], + Int(body_vortex_panels.totnode[]), + Int(body_vortex_panels.totpanel[]), + Int.(body_vortex_panels.nnode), + Int.(body_vortex_panels.npanel), body_vortex_panels.tangent, body_vortex_panels.controlpoint, - body_vortex_panels.endpanelidxs, + Int.(body_vortex_panels.endpanelidxs), idmaps.wake_node_ids_along_centerbody_wake_interface, idmaps.wake_node_ids_along_casing_wake_interface, idmaps.centerbody_panel_ids_along_centerbody_wake_interface, idmaps.duct_panel_ids_along_centerbody_wake_interface, + zpts, ) # extract surface velocity outputs @@ -320,17 +293,15 @@ function post_process( # Splits: vtan_casing_in, # tangent velocity along casing inside of duct body vtan_casing_out, # tangent velocity along casing outside of duct body - casing_zpts, # axial locations of control points along casing vtan_nacelle_in, # tangent velocity along nacelle inside of duct body vtan_nacelle_out, # tangent velocity along nacelle outside of duct body - nacelle_zpts, # axial locations of control points along nacelle vtan_centerbody_in, # tangent velocity inside of centerbody vtan_centerbody_out, # tangent velocity outside of centerbody - centerbody_zpts, # axial locations of control points along centerbody ) = vtan_tuple # - Surface Pressure on Bodies - # - cp_tuple = get_body_cps( + get_body_cps!( + cp_tuple, Vtan_in, Vtan_out, Gamr, @@ -344,6 +315,7 @@ function post_process( idmaps.id_of_first_centerbody_panel_aft_of_each_rotor, body_vortex_panels.controlpoint, body_vortex_panels.endpanelidxs, + zpts, ) # extract surface pressure outputs @@ -359,14 +331,20 @@ function post_process( ) = cp_tuple # - Calculate Thrust from Bodies - # - body_thrust, body_force_coeff = forces_from_pressure( - cp_in, cp_out, body_vortex_panels; rhoinf=rhoinf[1], Vref=Vref[1] + forces_from_pressure!( + body_thrust, + body_force_coefficient, + cp_in, + cp_out, + body_vortex_panels; + rhoinf=rhoinf[1], + Vref=Vref[1], ) # add thrust from trailing edge panels on bodies forces_from_TEpanels!( body_thrust, - body_force_coeff, + body_force_coefficient, cp_in, cp_out, body_vortex_panels; @@ -377,33 +355,41 @@ function post_process( ### --- TOTAL OUTPUTS --- ### # - Total Thrust - # - total_thrust = sum([rotor_inviscid_thrust'; rotor_viscous_thrust'; body_thrust]) + total_thrust[] = sum([rotor_inviscid_thrust'; rotor_viscous_thrust'; body_thrust]) # - Total Torque - # - total_torque = sum([rotor_inviscid_torque; rotor_viscous_torque]) + total_torque[] = sum([rotor_inviscid_torque; rotor_viscous_torque]) # - Total Power - # - total_power = sum([rotor_inviscid_power; rotor_viscous_power]) + total_power[] = sum([rotor_inviscid_power; rotor_viscous_power]) # - Total Efficiency - # - total_efficiency = get_total_efficiency(total_thrust, total_power, Vinf[1]) - ideal_efficiency = get_ideal_efficiency(total_thrust, rhoinf[1], Vinf[1], Rref[1]) + get_total_efficiency!(total_efficiency, total_thrust, total_power, Vinf[]) + ideal_efficiency[] = get_ideal_efficiency(total_thrust[], rhoinf[], Vinf[], Rref[]) # - Body Induced Efficiency - # - induced_efficiency = [ - get_induced_efficiency( - rotor_inviscid_thrust[ir], body_thrust, rotor_inviscid_power[ir], Vinf[1] - ) for ir in 1:nrotor - ] + get_induced_efficiency!( + induced_efficiency, + rotor_inviscid_thrust, + sum(body_thrust), + rotor_inviscid_power, + Vinf[], + ) # - Total Thrust and Torque Coefficients - # - total_CT, total_CQ, total_CP = tqpcoeff( - total_thrust, total_torque, total_power, rhoinf[1], Omega, Rref[1] + tqpcoeff!( + total_CT, + total_CQ, + total_CP, + total_thrust, + total_torque, + total_power, + rhoinf[1], + Omega, + Rref[1], ) outs = (; - # - States - # - # TODO: instead of states, call these something else so that the outputs are generalized # - Wake Values - # wake=(; panel_strengths=gamw), # - Body Values - # @@ -419,13 +405,13 @@ function post_process( cp_out, cp_casing_in, cp_casing_out, - casing_zpts, + zpts.casing_zpts, cp_nacelle_in, cp_nacelle_out, - nacelle_zpts, + zpts.nacelle_zpts, cp_centerbody_in, cp_centerbody_out, - centerbody_zpts, + zpts.centerbody_zpts, #individual body velocity contributions Vtot_in, Vtot_out, @@ -507,27 +493,233 @@ function post_process( ), ) - if write_outputs - write_data( - outs, - outfile; - output_tuple_name=output_tuple_name, - checkoutfileexists=checkoutfileexists, - ) - end - - return outs + if write_outputs + write_data( + outs, + outfile; + output_tuple_name=output_tuple_name, + checkoutfileexists=checkoutfileexists, + ) + end + + return outs +end + +###################################################################### +# # +# RESIDUAL FUNCTIONS # +# # +###################################################################### + +""" +""" +function run_residual!( + solver_options::TS, + converged_states, + state_dims, + solve_container_cache, + solve_container_cache_dims, + operating_point, + ivr, + ivw, + linsys, + blade_elements, + wakeK, + idmaps, +) where {TS<:ExternalSolverOptions} + + #= + NOTE: we want to get all the intermediate values available to user if desired. + The solve_containers cache will contain all the intermediate values after running the estimate states function. + =# + # - Separate out the state variables - # + vz_rotor, vtheta_rotor, Cm_wake = extract_state_variables( + solver_options, converged_states, state_dims + ) + + # - Extract and Reset Cache - # + # get cache vector of correct types + solve_container_cache_vec = @views PreallocationTools.get_tmp( + solve_container_cache, converged_states + ) + solve_containers = withdraw_solve_container_cache( + solver_options, solve_container_cache_vec, solve_container_cache_dims + ) + reset_containers!(solve_containers) #note: also zeros out state estimates + + # - Estimate New States - # + # currently has 280 allocations + estimate_states!( + solve_containers, + vz_rotor, + vtheta_rotor, + Cm_wake, + operating_point, + ivr, + ivw, + linsys, + blade_elements, + wakeK, + idmaps, + ) + + return (; vz_rotor, vtheta_rotor, Cm_wake, solve_containers...) +end + +""" +""" +function run_residual!( + solver_options::CSORSolverOptions, + converged_states, + state_dims, + solve_container_cache, + solve_container_cache_dims, + operating_point, + ivr, + ivw, + linsys, + blade_elements, + wakeK, + idmaps, +) + + #= + NOTE: we want to get all the intermediate values available to user if desired. + The solve_containers cache will contain all the intermediate values after running the insides of the residual function + =# + # - Separate out the state variables - # + Gamr, sigr, gamw = extract_state_variables(solver_options, converged_states, state_dims) + + # - Extract and Reset Cache - # + # get cache vector of correct types + solve_container_cache_vec = @views PreallocationTools.get_tmp( + solve_container_cache, converged_states + ) + solve_container_cache_vec .= 0 + solve_containers = withdraw_solve_container_cache( + solver_options, solve_container_cache_vec, solve_container_cache_dims + ) + + # - Run Residual - # + compute_CSOR_residual!( + zeros(2), + solver_options, + solve_containers, + Gamr, + sigr, + gamw, + operating_point, + ivr, + ivw, + linsys, + blade_elements, + wakeK, + idmaps; + verbose=false, + ) + + return (; Gamr, sigr, gamw, solve_containers...) +end + +###################################################################### +# # +# Velocity Functions # +# # +###################################################################### + +function get_body_tangential_velocities( + gamb, + gamw, + sigr, + ivb, + Vinf, + totnode, + totpanel, + nnode, + npanel, + tangent, + controlpoints, + endpanelidxs, + wake_panel_ids_along_centerbody_wake_interface, + wake_panel_ids_along_casing_wake_interface, + centerbody_panel_ids_along_centerbody_wake_interface, + duct_panel_ids_along_centerbody_wake_interface, + num_casing_panels, +) + TF = promote_type(eltype(gamb), eltype(gamw), eltype(sigr)) + # - initialize total velocity - # + Vtot = zeros(TF, 2, totpanel) + Vtot_in = similar(Vtot) .= 0 + Vtot_out = similar(Vtot) .= 0 + Vtan_in = similar(Vtot) .= 0 + Vtan_out = similar(Vtot) .= 0 + Vtot_prejump = similar(Vtot) .= 0 + vtot_body = similar(Vtot) .= 0 + duct_jump = similar(Vtot, (npanel[1],)) + centerbody_jump = similar(Vtot, (npanel[2],)) + body_jump_term = similar(Vtot) .= 0.0 + vtot_jump = similar(Vtot) .= 0.0 + vtot_wake = similar(Vtot) .= 0 + vtot_rotors = similar(Vtot) .= 0.0 + casing_zpts = zeros(TF, num_casing_panels) + vtan_casing_in = similar(casing_zpts) .= 0 + vtan_casing_out = similar(casing_zpts) .= 0 + nacelle_zpts = zeros(TF, npanels[1] - num_casing_panels) + vtan_nacelle_in = similar(nacelle_zpts) .= 0 + vtan_nacelle_out = similar(nacelle_zpts) .= 0 + centerbody_zpts = zeros(TF, npanels[2]) + vtan_centerbody_in = similar(centerbody_zpts) .= 0 + vtan_centerbody_out = similar(centerbody_zpts) .= 0 + + vtan_tuple = (; + # Totals and Components: + Vtan_in, + Vtot_in, + Vtan_out, + Vtot_out, + Vtot_prejump, + vtot_body, + duct_jump, + centerbody_jump, + body_jump_term, + vtot_jump, + vtot_wake, + vtot_rotors, + # Splits: + vtan_casing_in, + vtan_casing_out, + vtan_nacelle_in, + vtan_nacelle_out, + vtan_centerbody_in, + vtan_centerbody_out, + ) + + return get_body_tangential_velocities!( + vtan_tuple, + gamb, + gamw, + sigr, + ivb, + Vinf, + totnode, + totpanel, + nnode, + npanel, + tangent, + controlpoints, + endpanelidxs, + wake_panel_ids_along_centerbody_wake_interface, + wake_panel_ids_along_casing_wake_interface, + centerbody_panel_ids_along_centerbody_wake_interface, + duct_panel_ids_along_centerbody_wake_interface, + (; casing_zpts, nacelle_zpts, centerbody_zpts), + ) end -###################################################################### -# # -# Velocity Functions # -# # -###################################################################### - """ """ -function get_body_tangential_velocities( +function get_body_tangential_velocities!( + vtan_tuple, gamb, gamw, sigr, @@ -544,11 +736,11 @@ function get_body_tangential_velocities( wake_panel_ids_along_casing_wake_interface, centerbody_panel_ids_along_centerbody_wake_interface, duct_panel_ids_along_centerbody_wake_interface, + zpts, ) # - setup - # nws, nrotor = size(sigr) - TF = promote_type(eltype(gamb), eltype(gamw), eltype(sigr)) (; v_bb, v_br, v_bw) = ivb # rename for convenience @@ -557,99 +749,117 @@ function get_body_tangential_velocities( whi = centerbody_panel_ids_along_centerbody_wake_interface wdi = duct_panel_ids_along_centerbody_wake_interface - # TODO also consider including the body wakes here as well. + (; + # Totals and Components: + Vtan_in, + Vtot_in, + Vtan_out, + Vtot_out, + Vtot_prejump, + vtot_body, + duct_jump, + centerbody_jump, + body_jump_term, + vtot_jump, + vtot_wake, + vtot_rotors, + # Splits: + vtan_casing_in, + vtan_casing_out, + vtan_nacelle_in, + vtan_nacelle_out, + vtan_centerbody_in, + vtan_centerbody_out, + ) = vtan_tuple - # - initialize total velocity - # - Vtot = zeros(TF, 2, totpanel) + (; casing_zpts, nacelle_zpts, centerbody_zpts) = zpts - # - Velocity Contributions from body - # + # TODO also consider including the body wakes here as well. - vtot_body = similar(Vtot) .= 0 + # - Velocity Contributions from body - # for (i, vt) in enumerate(eachrow(vtot_body)) vt .= @views v_bb[:, :, i] * gamb[1:size(v_bb, 2)] end - Vtot .+= vtot_body + Vtot_out .+= vtot_body # - Velocity Contributions from wake - # - vtot_wake = similar(Vtot) .= 0 for (i, vt) in enumerate(eachrow(vtot_wake)) vt .= @views v_bw[:, :, i] * gamw end - Vtot .+= vtot_wake # opposite sign from linear solve + Vtot_out .+= vtot_wake # opposite sign from linear solve # - Velocity Contributions from rotors - # - vtot_rotors = similar(Vtot) .= 0.0 for jrotor in 1:nrotor rotorrange = (nws * (jrotor - 1) + 1):(nws * jrotor) for (i, vt) in enumerate(eachrow(vtot_rotors)) vt .+= @views v_br[:, rotorrange, i] * sigr[rotorrange] end end - Vtot .+= vtot_rotors # opposite sign from linear solve + Vtot_out .+= vtot_rotors # opposite sign from linear solve # - Influence from Freestream - # - Vtot[1, :] .+= Vinf # opposite sign from linear solve - Vtot_prejump = copy(Vtot) + Vtot_out[1, :] .+= Vinf # opposite sign from linear solve + Vtot_prejump .= copy(Vtot_out) # - Add in Jump Term - # # duct - jumpduct = (gamb[1:(npanel[1])] + gamb[2:(nnode[1])]) / 2 + duct_jump .= @views (gamb[1:(npanel[1])] + gamb[2:(nnode[1])]) / 2 # wake panels interfacing with duct - jumpduct[wdi] .+= (gamw[dwi[1]:(dwi[end] - 1)] + gamw[(dwi[1] + 1):dwi[end]]) / 2.0 + duct_jump[wdi] .+= @views (gamw[dwi[1]:(dwi[end] - 1)] + gamw[(dwi[1] + 1):dwi[end]]) / + 2.0 # center body panels - jumphub = (gamb[(nnode[1] + 1):(totnode - 1)] + gamb[(nnode[1] + 2):(totnode)]) / 2.0 + centerbody_jump .= @views ( + gamb[(nnode[1] + 1):(totnode - 1)] + gamb[(nnode[1] + 2):(totnode)] + ) / 2.0 # wake panels interfacing with center body - jumphub[whi] .+= (gamw[hwi[1]:(hwi[end] - 1)] + gamw[(hwi[1] + 1):hwi[end]]) / 2.0 + centerbody_jump[whi] .+= @views ( + gamw[hwi[1]:(hwi[end] - 1)] + gamw[(hwi[1] + 1):hwi[end]] + ) / 2.0 - body_jump_term = [jumpduct; jumphub] + body_jump_term[1:length(duct_jump)] .= duct_jump + body_jump_term[(length(duct_jump) + 1):end] .= centerbody_jump - vtot_jump = similar(Vtot) .= 0.0 for (vt, tan) in zip(eachrow(vtot_jump), eachrow(tangent)) vt .+= body_jump_term .* tan ./ 2.0 end # assign velocities to each side of the panel - Vtot_out = Vtot .- vtot_jump # outer side of boundary - Vtot_in = Vtot .+ vtot_jump # inner side of boundary + Vtot_in .= Vtot_out .+ vtot_jump # inner side of boundary + Vtot_out .-= vtot_jump # outer side of boundary # Get the magnitude of the sum of the velocities and this is the surface velocity since the body velocities have been solved to eliminate the normal components in the summed velocities - Vtan_out = sqrt.(Vtot_out[1, :] .^ 2 .+ Vtot_out[2, :] .^ 2) - Vtan_in = sqrt.(Vtot_in[1, :] .^ 2 .+ Vtot_in[2, :] .^ 2) + Vtan_out .= @views sqrt.(Vtot_out[1, :] .^ 2 .+ Vtot_out[2, :] .^ 2) + Vtan_in .= @views sqrt.(Vtot_in[1, :] .^ 2 .+ Vtot_in[2, :] .^ 2) # - Split Velocities associates with inner and outer duct and hub - # # total tangential velocities - vtan_casing_out, vtan_nacelle_out, vtan_centerbody_out, casing_zpts, nacelle_zpts, centerbody_zpts = split_bodies( - Vtan_out, controlpoints, endpanelidxs - ) - vtan_casing_in, vtan_nacelle_in, vtan_centerbody_in, casing_zpts, nacelle_zpts, centerbody_zpts = split_bodies( - Vtan_in, controlpoints, endpanelidxs - ) - - return (; - # Totals and Components: - Vtan_in, - Vtot_in, - Vtan_out, - Vtot_out, - Vtot_prejump, - vtot_body, - vtot_jump, - vtot_wake, - vtot_rotors, - # Splits: - vtan_casing_in, + split_bodies!( vtan_casing_out, - casing_zpts, - vtan_nacelle_in, vtan_nacelle_out, + vtan_centerbody_out, + casing_zpts, nacelle_zpts, + centerbody_zpts, + Vtan_out, + controlpoints, + endpanelidxs, + ) + split_bodies!( + vtan_casing_in, + vtan_nacelle_in, vtan_centerbody_in, - vtan_centerbody_out, + casing_zpts, + nacelle_zpts, centerbody_zpts, + Vtan_in, + controlpoints, + endpanelidxs, ) + + return vtan_tuple end """ @@ -713,9 +923,14 @@ end Calculate steady pressure coefficient """ function steady_cp(Vs, Vinf, Vref) - return (Vinf^2 .- Vs .^ 2) / Vref^2 + cp = similar(Vs) .= 0 + return steady_cp!(cp, Vs, Vinf, Vref) end +function steady_cp!(cp, Vs, Vinf, Vref) + cp .= (Vinf^2 .- Vs .^ 2) / Vref^2 + return cp +end """ only used in post-process for cp. expression not in dfdc theory, comes from source code. @@ -818,11 +1033,6 @@ function calculate_bodywake_delta_cp(Gamr, sigr, Cz_rotor, Vref, Omega, B, r; bo return deltacp end -""" -calculate pressure coefficient distribution on duct/hub walls -formulation taken from DFDC source code. TODO: derive where the expressions came from. - -""" function get_body_cps( Vtan_in, Vtan_out, @@ -837,11 +1047,82 @@ function get_body_cps( hidr, controlpoints, endpanelidxs, + zpts, +) + cp_in = similar(Vtan_in) .= 0 + cp_out = similar(Vtan_in) .= 0 + cp_casing_in = similar(Vtan_in, size(zpts.casing_zpts)) .= 0 + cp_casing_out = similar(Vtan_in, size(zpts.casing_zpts)) .= 0 + cp_nacelle_in = similar(Vtan_in, size(zpts.nacelle_zpts)) .= 0 + cp_nacelle_out = similar(Vtan_in, size(zpts.nacelle_zpts)) .= 0 + cp_centerbody_in = similar(Vtan_in, size(zpts.centerbody_zpts)) .= 0 + cp_centerbody_out = similar(Vtan_in, size(zpts.centerbody_zpts)) .= 0 + + cp_tuple = (; + cp_in, + cp_out, + cp_casing_in, + cp_casing_out, + cp_nacelle_in, + cp_nacelle_out, + cp_centerbody_in, + cp_centerbody_out, + ) + + return get_body_cps!( + cp_tuple, + Vtan_in, + Vtan_out, + Gamr, + sigr, + Cz_rotor, + Vinf, + Vref, + B, + Omega, + didr, + hidr, + controlpoints, + endpanelidxs, + zpts, + ) +end + +function get_body_cps!( + cp_tuple, + Vtan_in, + Vtan_out, + Gamr, + sigr, + Cz_rotor, + Vinf, + Vref, + B, + Omega, + didr, + hidr, + controlpoints, + endpanelidxs, + zpts, ) + # rename for convenience + (; + cp_in, # surface pressure along inside of bodies + cp_out, # surface pressure along outside of bodies + cp_casing_in, # surface pressure along inside of casing + cp_nacelle_in, # surface pressure along inside of nacell + cp_centerbody_in, # surface pressure along inside of centerbody + cp_casing_out, # surface pressure along outside of casing + cp_nacelle_out, # surface pressure along outside of nacelle + cp_centerbody_out, # surface pressure along outside of centerbody + ) = cp_tuple + + (; casing_zpts, nacelle_zpts, centerbody_zpts) = zpts + # - Calculate standard pressure coefficient expression - # - cp_in = steady_cp(Vtan_in, Vinf, Vref) - cp_out = steady_cp(Vtan_out, Vinf, Vref) + steady_cp!(cp_in, Vtan_in, Vinf, Vref) + steady_cp!(cp_out, Vtan_out, Vinf, Vref) # - add the change in Cp on the walls due to enthalpy, entropy, and vtheta - # calculate_body_delta_cp!( @@ -849,26 +1130,30 @@ function get_body_cps( ) # - Split body strengths into inner/outer duct and hub - # - cp_casing_in, cp_nacelle_in, cp_centerbody_in, casing_zpts, nacelle_zpts, centerbody_zpts = split_bodies( - cp_in, controlpoints, endpanelidxs - ) - cp_casing_out, cp_nacelle_out, cp_centerbody_out, _, _, _ = split_bodies( - cp_out, controlpoints, endpanelidxs - ) - - return (; - cp_in, - cp_out, + split_bodies!( cp_casing_in, cp_nacelle_in, cp_centerbody_in, casing_zpts, nacelle_zpts, centerbody_zpts, + cp_in, + controlpoints, + endpanelidxs, + ) + split_bodies!( cp_casing_out, cp_nacelle_out, cp_centerbody_out, + casing_zpts, + nacelle_zpts, + centerbody_zpts, + cp_out, + controlpoints, + endpanelidxs, ) + + return cp_tuple end """ @@ -920,31 +1205,36 @@ end Calculate dimensional and non-dimensional axial force on a single body """ function forces_from_pressure(cp_in, cp_out, panels; rhoinf=1.225, Vref=1.0) + # - initialize - # + cfx = zeros(eltype(cp_out), Int(panels.nbodies[])) # axial force coefficient (all others are zero for axisymmetric case) + CFx = similar(cfx) .= 0 + + return forces_from_pressure!(CFx, cfx, cp_in, cp_out, panels; rhoinf=rhoinf, Vref=Vref) +end + +function forces_from_pressure!(CFx, cfx, cp_in, cp_out, panels; rhoinf=1.225, Vref=1.0) # - rename for convenience - # #just want x-component of normals since it's axisymmetric - ns = panels.normal[1, :] + ns = @view(panels.normal[1, :]) #radial positions - rs = panels.controlpoint[2, :] + rs = @view(panels.controlpoint[2, :]) #panel lengths ds = panels.influence_length - # - initialize - # - cfx = zeros(eltype(cp_out), panels.nbodies) # axial force coefficient (all others are zero for axisymmetric case) - # for each body - for ib in 1:(panels.nbodies) + for ib in 1:(Int(panels.nbodies[])) # - rectangular integration due to constant panel strengths. - # - for ip in panels.endpanelidxs[1, ib]:panels.endpanelidxs[2, ib] + for ip in Int.(panels.endpanelidxs[1, ib]:panels.endpanelidxs[2, ib]) cfx[ib] += (cp_out[ip] - cp_in[ip]) * ns[ip] * ds[ip] * 2.0 * pi * rs[ip] end end #dimensionalize - q = 0.5 * rhoinf * Vref^2 + CFx .= cfx .* 0.5 * rhoinf * Vref^2 #note, thrust is in negative x-direction - return cfx .* q, cfx + return CFx, cfx end """ @@ -957,18 +1247,23 @@ function forces_from_TEpanels!( #dimensionalize q = 0.5 * rhoinf * Vref^2 - for i in 1:(panels.nbodies) + for i in 1:(Int(panels.nbodies[])) if panels.tenode[i, 2, 2] <= eps() # if it's the hub, don't average the first and last, but rather just the last - cpi = cp_in[panels.endpanelidxs[i, 2]] - cpo = cp_out[panels.endpanelidxs[i, 2]] + cpi = cp_in[Int(panels.endpanelidxs[i, 2])] + cpo = cp_out[Int(panels.endpanelidxs[i, 2])] else # if it's the duct, then average the first and last panel cpi = - 0.5 * (cp_in[panels.endpanelidxs[1, i]] + cp_in[panels.endpanelidxs[2, i]]) + 0.5 * ( + cp_in[Int(panels.endpanelidxs[1, i])] + + cp_in[Int(panels.endpanelidxs[2, i])] + ) cpo = - 0.5 * - (cp_out[panels.endpanelidxs[1, i]] + cp_out[panels.endpanelidxs[2, i]]) + 0.5 * ( + cp_out[Int(panels.endpanelidxs[1, i])] + + cp_out[Int(panels.endpanelidxs[2, i])] + ) end r = 0.5 * sum(panels.tenode[i, :, 2]) @@ -1001,12 +1296,21 @@ end ###################################################################### function inviscid_rotor_thrust(Ctheta_rotor, Gamma_tilde, rotor_panel_length, rhoinf) - - # problem dimensions - nr, nrotor = size(Gamma_tilde) - # initialize dTi = similar(Gamma_tilde) .= 0.0 + Tinv = zeros(eltype(Gamma_tilde), size(Gamma_tilde, 2)) + + return inviscid_rotor_thrust!( + Tinv, dTi, Ctheta_rotor, Gamma_tilde, rotor_panel_length, rhoinf + ) +end + +function inviscid_rotor_thrust!( + Tinv, dTi, Ctheta_rotor, Gamma_tilde, rotor_panel_length, rhoinf +) + + # problem dimensions + nr, nrotor = size(dTi) for irotor in 1:nrotor for ir in 1:nr @@ -1020,7 +1324,7 @@ function inviscid_rotor_thrust(Ctheta_rotor, Gamma_tilde, rotor_panel_length, rh end #sum the section thrust - Tinv = sum(dTi; dims=1) + Tinv .= sum(dTi; dims=1) return Tinv, dTi end @@ -1028,22 +1332,35 @@ end function viscous_rotor_thrust( Cz_rotor, Cmag_rotor, B, chord, rotor_panel_length, cd, rhoinf ) - - # get dimensions - nr, nrotor = size(Cz_rotor) - #initialize dTv = similar(Cz_rotor) .= 0.0 + Tvisc = zeros(eltype(Cz_rotor), size(Cz_rotor, 2)) + + return viscous_rotor_thrust!( + Tvisc, dTv, Cz_rotor, Cmag_rotor, B, chord, rotor_panel_length, cd, rhoinf + ) +end + +function viscous_rotor_thrust!( + Tvisc, dTv, Cz_rotor, Cmag_rotor, B, chord, rotor_panel_length, cd, rhoinf +) + + # get dimensions + nr, nrotor = size(dTv) for irotor in 1:nrotor for ir in 1:nr - hrwc = 0.5 * rhoinf * Cmag_rotor[ir, irotor] * chord[ir, irotor] - bdr = B[irotor] * rotor_panel_length[ir, irotor] - dTv[ir, irotor] = -hrwc * cd[ir, irotor] * Cz_rotor[ir, irotor] * bdr + # hrwc = 0.5 * rhoinf * Cmag_rotor[ir, irotor] * chord[ir, irotor] + # bdr = B[irotor] * rotor_panel_length[ir, irotor] + dTv[ir, irotor] = + -(0.5 * rhoinf * Cmag_rotor[ir, irotor] * chord[ir, irotor]) * + cd[ir, irotor] * + Cz_rotor[ir, irotor] * + (B[irotor] * rotor_panel_length[ir, irotor]) end end - Tvisc = sum(dTv; dims=1) + Tvisc .= sum(dTv; dims=1) return Tvisc, dTv end @@ -1051,21 +1368,34 @@ end function inviscid_rotor_torque( Cz_rotor, rotor_panel_center, rotor_panel_length, Gamma_tilde, rhoinf ) - - # dimensions - nr, nrotor = size(Gamma_tilde) - # initialize dQi = similar(Gamma_tilde) .= 0.0 + Qinv = zeros(eltype(Gamma_tilde), size(Gamma_tilde, 2)) + + return inviscid_rotor_torque!( + Qinv, dQi, Cz_rotor, rotor_panel_center, rotor_panel_length, Gamma_tilde, rhoinf + ) +end + +function inviscid_rotor_torque!( + Qinv, dQi, Cz_rotor, rotor_panel_center, rotor_panel_length, Gamma_tilde, rhoinf +) + + # dimensions + nr, nrotor = size(dQi) for irotor in 1:nrotor for ir in 1:nr - rdr = rotor_panel_center[ir, irotor] * rotor_panel_length[ir, irotor] - dQi[ir, irotor] = rhoinf * Gamma_tilde[ir, irotor] * Cz_rotor[ir, irotor] * rdr + # rdr = rotor_panel_center[ir, irotor] * rotor_panel_length[ir, irotor] + dQi[ir, irotor] = + rhoinf * + Gamma_tilde[ir, irotor] * + Cz_rotor[ir, irotor] * + (rotor_panel_center[ir, irotor] * rotor_panel_length[ir, irotor]) end end - Qinv = sum(dQi; dims=1) + Qinv .= sum(dQi; dims=1) return Qinv, dQi end @@ -1081,11 +1411,6 @@ function viscous_rotor_torque( rhoinf; TF=nothing, ) - - # dimensions - nr, nrotor = size(Wtheta_rotor) - - # initialize if isnothing(TF) TF = promote_type( eltype(Wtheta_rotor), @@ -1096,28 +1421,84 @@ function viscous_rotor_torque( ) end - dQv = zeros(TF, size(chord)) + dQv = zeros(TF, size(Cmag_rotor)) + Qvisc = zeros(TF, size(B)) + + return viscous_rotor_torque!( + Qvisc, + dQv, + Wtheta_rotor, + Cmag_rotor, + B, + chord, + rotor_panel_center, + rotor_panel_length, + cd, + rhoinf; + TF=nothing, + ) +end + +function viscous_rotor_torque!( + Qvisc, + dQv, + Wtheta_rotor, + Cmag_rotor, + B, + chord, + rotor_panel_center, + rotor_panel_length, + cd, + rhoinf; + TF=nothing, +) + + # dimensions + nr, nrotor = size(dQv) + + # initialize for irotor in 1:nrotor for ir in 1:nr - hrwc = 0.5 * rhoinf * Cmag_rotor[ir, irotor] * chord[ir, irotor] - brdr = - B[irotor] * rotor_panel_center[ir, irotor] * rotor_panel_length[ir, irotor] - dQv[ir, irotor] = -hrwc * cd[ir, irotor] * Wtheta_rotor[ir, irotor] * brdr + # hrwc = 0.5 * rhoinf * Cmag_rotor[ir, irotor] * chord[ir, irotor] + # brdr = + # B[irotor] * rotor_panel_center[ir, irotor] * rotor_panel_length[ir, irotor] + dQv[ir, irotor] = + -(0.5 * rhoinf * Cmag_rotor[ir, irotor] * chord[ir, irotor]) * + cd[ir, irotor] * + Wtheta_rotor[ir, irotor] * + ( + B[irotor] * + rotor_panel_center[ir, irotor] * + rotor_panel_length[ir, irotor] + ) end end - Qvisc = sum(dQv; dims=1) + Qvisc .= sum(dQv; dims=1) return Qvisc, dQv end -function inviscid_rotor_power(Qinv, Omega) - return Qinv .* Omega +function rotor_power(Q, dQ, Omega) + dP = similar(dQ) .= 0.0 + P = similar(Q) .= 0.0 + + return rotor_power!(P, dP, Q, dQ, Omega) end -function viscous_rotor_power(Qvisc, Omega) - return Qvisc .* Omega +function rotor_power!(P, dP, Q, dQ, Omega) + nr, nrotor = size(dP) + + for irotor in 1:nrotor + for ir in 1:nr + dP[ir, irotor] = dQ[ir, irotor] * Omega[irotor] + end + + P[irotor] = Q[irotor] .* Omega[irotor] + end + + return P, dP end function get_total_efficiency(total_thrust, total_power, Vinf) @@ -1125,6 +1506,10 @@ function get_total_efficiency(total_thrust, total_power, Vinf) eta = zeros(TF, length(total_thrust)) + return get_total_efficiency!(eta, total_thrust, total_power, Vinf) +end + +function get_total_efficiency!(eta, total_thrust, total_power, Vinf) for i in 1:length(total_thrust) if Vinf <= 0.0 || total_power[i] < eps() || total_thrust[i] <= 0.0 #do nothing, efficiency can't physically be negative or infinite. @@ -1137,17 +1522,29 @@ function get_total_efficiency(total_thrust, total_power, Vinf) end function get_induced_efficiency(Tinv, Tduct, Pinv, Vinf) - if Vinf <= 0.0 || Pinv <= 0.0 - return 0.0 - else - return Vinf * (Tinv .+ Tduct) ./ Pinv + TF = promote_type(eltype(Tinv), eltype(Pinv), eltpye(Tduct)) + eta_inv = zeros(TF, size(Tinv)) + return get_induced_efficiency!(eta_inv, Tinv, Tduct, Pinv, Vinf) +end + +function get_induced_efficiency!(eta_inv, Tinv, Tduct, Pinv, Vinf) + for (e, ti, p) in zip(eachrow(eta_inv), Tinv, Pinv) + if Vinf <= 0.0 || p <= 0.0 + e[1] = 0.0 + else + e[1] = Vinf * (ti + Tduct) / p + end end + return eta_inv end function get_ideal_efficiency(total_thrust, rhoinf, Vinf, Rref) if Vinf != 0.0 - TC = total_thrust / (0.5 * rhoinf * Vinf^2 * pi * Rref^2) - return 2.0 / (1.0 + sqrt(max(TC, -1.0) + 1.0)) + # TC = total_thrust / (0.5 * rhoinf * Vinf^2 * pi * Rref^2) + return 2.0 / ( + 1.0 + + sqrt(max((total_thrust / (0.5 * rhoinf * Vinf^2 * pi * Rref^2)), -1.0) + 1.0) + ) else return 0.0 end @@ -1158,25 +1555,28 @@ function tqpcoeff(thrust, torque, power, rhoinf, Omega, Rref) CT = zeros(T, length(Omega)) CQ = zeros(T, length(Omega)) CP = zeros(T, length(Omega)) + return tqpcoeff!(CT, CQ, CP, thrust, torque, power, rhoinf, Omega, Rref) +end +function tqpcoeff!(CT, CQ, CP, thrust, torque, power, rhoinf, Omega, Rref) for (i, o) in enumerate(Omega) if isapprox(o, 0.0) CT[i] = CQ[i] = CP[i] = 0.0 else # reference diameter - D = 2.0 * Rref + # D = 2.0 * Rref # rototion in rev per second - n = o / (2.0 * pi) + # n = o / (2.0 * pi) # thrust coefficient - CT[i] = thrust[i] / (rhoinf * n^2 * D^4) + CT[i] = thrust[i] / (rhoinf * (o / (2.0 * pi))^2 * (2.0 * Rref)^4) # torque coefficient - CQ[i] = torque[i] / (rhoinf * n^2 * D^5) + CQ[i] = torque[i] / (rhoinf * (o / (2.0 * pi))^2 * (2.0 * Rref)^5) # power coefficient - CP[i] = power[i] / (rhoinf * n^3 * D^5) + CP[i] = power[i] / (rhoinf * (o / (2.0 * pi))^3 * (2.0 * Rref)^5) end end @@ -1184,37 +1584,50 @@ function tqpcoeff(thrust, torque, power, rhoinf, Omega, Rref) end function get_blade_loads(Wmag_rotor, beta1, cl, cd, chords, rhoinf)#, Rhub, Rtip, rotor_panel_centers,rotor_panel_lengths ,B, Omega) - - # dimensions - nr, nrotor = size(Wmag_rotor) - # initialize Np = similar(Wmag_rotor) .= 0.0 Tp = similar(Wmag_rotor) .= 0.0 + return get_blade_loads!(Np, Tp, Wmag_rotor, beta1, cl, cd, chords, rhoinf)#, Rhub, Rtip, rotor_panel_centers,rotor_panel_lengths ,B, Omega) +end + +function get_blade_loads!(Np, Tp, Wmag_rotor, beta1, cl, cd, chords, rhoinf, cache)#, Rhub, Rtip, rotor_panel_centers,rotor_panel_lengths ,B, Omega) + + # dimensions + nr, nrotor = size(Np) for irotor in 1:nrotor for ir in 1:nr # rename for convenience - cphi = cos(beta1[ir, irotor]) - sphi = sin(beta1[ir, irotor]) + cache.cphi[ir, irotor] = cos(beta1[ir, irotor]) + cache.sphi[ir, irotor] = sin(beta1[ir, irotor]) # resolve lift and drag into normal and tangential coefficients - cn = cl[ir, irotor] * cphi - cd[ir, irotor] * sphi - ct = cl[ir, irotor] * sphi + cd[ir, irotor] * cphi + cache.cn[ir, irotor] = + cl[ir, irotor] * cache.cphi[ir, irotor] - + cd[ir, irotor] * cache.sphi[ir, irotor] + cache.ct[ir, irotor] = + cl[ir, irotor] * cache.sphi[ir, irotor] + + cd[ir, irotor] * cache.cphi[ir, irotor] # get the normal and tangential loads per unit length N' and T' Np[ir, irotor] = - cn * 0.5 * rhoinf * Wmag_rotor[ir, irotor]^2 * chords[ir, irotor] + cache.cn[ir, irotor] * + 0.5 * + rhoinf * + Wmag_rotor[ir, irotor]^2 * + chords[ir, irotor] Tp[ir, irotor] = - ct * 0.5 * rhoinf * Wmag_rotor[ir, irotor]^2 * chords[ir, irotor] + cache.ct[ir, irotor] * + 0.5 * + rhoinf * + Wmag_rotor[ir, irotor]^2 * + chords[ir, irotor] end end # Npfull = [zeros(nrotor)'; Np; zeros(nrotor)'] # Tpfull = [zeros(nrotor)'; Tp; zeros(nrotor)'] - # TODO: consider comparing this with DFDC versions for them - ## -- Integrate Loads to get Thrust and Torque # add hub/tip for complete integration. loads go to zero at hub/tip. # rfull = [Rhub; rotor_panel_centers; Rtip] diff --git a/src/preprocess/geometry/panel_geometry.jl b/src/preprocess/geometry/panel_geometry.jl index c40bcee4..5a8d3953 100644 --- a/src/preprocess/geometry/panel_geometry.jl +++ b/src/preprocess/geometry/panel_geometry.jl @@ -6,11 +6,11 @@ """ avoid issues with single bodies being input not as vectors """ -function generate_panels(coordinates::Matrix{TF}; kwargs...) where {TF} +function generate_panels(coordinates::AbstractMatrix{TF}; kwargs...) where {TF} return generate_panels([coordinates]; kwargs...) end -function generate_panels!(panels, coordinates::Matrix{TF}; kwargs...) where {TF} +function generate_panels!(panels, coordinates::AbstractMatrix{TF}; kwargs...) where {TF} return generate_panels!(panels, [coordinates]; kwargs...) end @@ -33,54 +33,49 @@ function generate_panels( # number of panels to generate for each body npanel = nnode .- 1 # number of bodies - nbodies = length(npanel) + nbodies = [length(npanel)] # total number of panels in system - totpanel = sum(npanel) + totpanel = [sum(npanel)] # total number of nodes in system totnode = totpanel + nbodies # - Initialize Outputs - # # control points - controlpoint = zeros(TF, 2, totpanel) # x-r, panel + controlpoint = zeros(TF, 2, totpanel[]) # x-r, panel # nodes - node = zeros(TF, 2, totnode) # x-r, node + node = zeros(TF, 2, totnode[]) # x-r, node # node map - nodemap = zeros(Int, 2, totpanel) # node, x-r + nodemap = zeros(Int, 2, totpanel[]) # node, x-r # first and last nodes of bodies, rotor, or wakes #TODO: check where this is used and update order of indices if needed - endnodes = zeros(TF, nbodies, 2, 2) # TE, upper-lower, x-r + endnodes = zeros(TF, nbodies[], 2, 2) # TE, upper-lower, x-r # indices of endnodess - endnodeidxs = ones(Int, 2, nbodies) # lower idx, upper idx, lower or upper + endnodeidxs = ones(Int, 2, nbodies[]) # lower idx, upper idx, lower or upper # indices of endpanels - endpanelidxs = ones(Int, 2, nbodies) # lower idx, upper idx, lower or upper + endpanelidxs = ones(Int, 2, nbodies[]) # lower idx, upper idx, lower or upper # panel lengths - influence_length = zeros(TF, totpanel) + influence_length = zeros(TF, totpanel[]) # panel unit normals - normal = zeros(TF, 2, totpanel) + normal = zeros(TF, 2, totpanel[]) # panel unit tangents - tangent = zeros(TF, 2, totpanel) + tangent = zeros(TF, 2, totpanel[]) # internal control point - itcontrolpoint = zeros(TF, 2, nbodies) - itnormal = zeros(TF, 2, nbodies) + itcontrolpoint = zeros(TF, 2, nbodies[]) + itnormal = zeros(TF, 2, nbodies[]) #note: unused, but required input later - ittangent = zeros(TF, 2, nbodies) + ittangent = zeros(TF, 2, nbodies[]) # trailing edge gap - tenode = zeros(TF, nbodies, 2, 2) # body, node1-2, x-r - tenormal = zeros(TF, 2, nbodies) #body, x-r - teinfluence_length = zeros(TF, nbodies) + tenode = zeros(TF, nbodies[], 2, 2) # body, node1-2, x-r + tenormal = zeros(TF, 2, nbodies[]) #body, x-r + teinfluence_length = zeros(TF, nbodies[]) teadjnodeidxs = similar(endnodeidxs) .= 1 #can't be the same as endpoints, because we may have repeated values for non-duct bodies - tendotn = zeros(TF, 2, nbodies) #bodies, node1,2 - tencrossn = zeros(TF, 2, nbodies) #bodies, node1,2 + tendotn = zeros(TF, 2, nbodies[]) #bodies, node1,2 + tencrossn = zeros(TF, 2, nbodies[]) #bodies, node1,2 - # count number of prescribed nodes - npnid = 0 - for c in coordinates - npnid += count(x -> abs(x) <= eps(), c[2, :]) - end - prescribednodeidxs = zeros(Int, npnid) + prescribednodeidxs = zeros(Int, nbodies[]) panels = (; controlpoint, @@ -122,7 +117,7 @@ end function generate_panels!( panels, - coordinates::Vector{Matrix{TF}}; + coordinates::AbstractVector{TF}; itcpshift=0.05, axistol=1e-15, tegaptol=1e1 * eps(), @@ -158,6 +153,20 @@ function generate_panels!( totpanel, ) = panels + if iszero(nnode) + ## -- SETUP -- ## + # number of nodes according to coordinates + nnode .= [length(eachcol(c)) for c in coordinates] + # number of panels to generate for each body + npanel .= nnode .- 1 + # number of bodies + nbodies .= length(npanel) + # total number of panels in system + totpanel .= sum(npanel) + # total number of nodes in system + totnode .= totpanel + nbodies + end + # initialize index for entire array pidx = 1 # panel index nidx = 1 # node index (+1 extra after each body) @@ -174,11 +183,11 @@ function generate_panels!( @assert all(r -> r >= 0.0, r) "Some coordinates have negative radial components." ## -- Loop Through Coordinates -- ## - for ip in 1:npanel[ib] + for ip in 1:Int(npanel[ib]) # Get nodes (panel edges) - node[:, nidx] = [x[ip]; r[ip]] - node[:, nidx + 1] = [x[ip + 1]; r[ip + 1]] + node[:, nidx] .= [x[ip]; r[ip]] + node[:, nidx + 1] .= [x[ip + 1]; r[ip + 1]] nodemap[:, pidx] .= [nidx; nidx + 1] # Calculate control point (panel center) @@ -199,7 +208,7 @@ function generate_panels!( endnodes[ib, 1, :] = [x[ip] r[ip]] endnodeidxs[1, ib] = nidx endpanelidxs[1, ib] = pidx - elseif ip == npanel[ib] + elseif ip == Int(npanel[ib]) #upper endnodes[ib, 2, :] = [x[ip + 1] r[ip + 1]] endnodeidxs[2, ib] = nidx + 1 @@ -252,7 +261,12 @@ function generate_panels!( # - Prescribed Nodes - # # Save the node index for nodes that are on the axis and need to be prescribed. - prescribednodeidxs[:] = findall(x -> abs(x) <= eps(), node[2, :]) + if isbody + prescribednodeidxs[1] = findfirst(x -> abs(x) <= eps(), node[2, :]) + prescribednodeidxs[2] = node[2, end] <= eps() ? size(node, 2) : 0 + else + prescribednodeidxs .= 0 + end return panels end @@ -270,11 +284,11 @@ function def_it_panel!( ) #TODO: maybe move this into it's own in-place function - for ib in 1:nbodies + for ib in 1:Int(nbodies[]) #rename for convenience - p1id = endnodeidxs[1, ib] - pnid = endnodeidxs[2, ib] + p1id = Int(endnodeidxs[1, ib]) + pnid = Int(endnodeidxs[2, ib]) # get node coordinates n1 = node[:, p1id] #first node on panel 1 (of body ib) @@ -315,9 +329,10 @@ function def_te_panel!( endnodes, endnodeidxs, ) - for ib in 1:nbodies + for ib in 1:Int(nbodies[]) # check if signs of x-tangents are the same - if sign(tangent[1, endpanelidxs[1, ib]]) != sign(tangent[1, endpanelidxs[2, ib]]) + if sign(tangent[1, Int(endpanelidxs[1, ib])]) != + sign(tangent[1, Int(endpanelidxs[2, ib])]) # if not: it's a duct # set first node to last endnode, tenode[ib, 1, :] .= endnodes[ib, 2, :] @@ -326,13 +341,17 @@ function def_te_panel!( # set normal as above tenormal[:, ib] = get_panel_normal(get_r(tenode[ib, 1, :], tenode[ib, 2, :])...) # set node id's to the endnode ids - teadjnodeidxs[1, ib] = endnodeidxs[2, ib] - teadjnodeidxs[2, ib] = endnodeidxs[1, ib] + teadjnodeidxs[1, ib] = Int(endnodeidxs[2, ib]) + teadjnodeidxs[2, ib] = Int(endnodeidxs[1, ib]) # set dots and crosses for each adjacent panel - tendotn[1, ib] = dot(tenormal[:, ib], normal[:, endpanelidxs[2, ib]]) - tendotn[2, ib] = dot(tenormal[:, ib], normal[:, endpanelidxs[1, ib]]) - tencrossn[1, ib] = cross2mag(normal[:, endpanelidxs[2, ib]], tenormal[:, ib]) - tencrossn[2, ib] = cross2mag(normal[:, endpanelidxs[1, ib]], tenormal[:, ib]) + tendotn[1, ib] = dot(tenormal[:, ib], normal[:, Int(endpanelidxs[2, ib])]) + tendotn[2, ib] = dot(tenormal[:, ib], normal[:, Int(endpanelidxs[1, ib])]) + tencrossn[1, ib] = cross2mag( + normal[:, Int(endpanelidxs[2, ib])], tenormal[:, ib] + ) + tencrossn[2, ib] = cross2mag( + normal[:, Int(endpanelidxs[1, ib])], tenormal[:, ib] + ) else # if so: it's anything else (assuming the hub nose doesn't curve inward... @@ -343,13 +362,15 @@ function def_te_panel!( # set normal to parallel to axis tenormal[:, ib] = [1.0, 0.0] # set both adjacent node id's to the last endnode id - # teadjnodeidxs[ib, :] .= endnodeidxs[2, ib] - teadjnodeidxs[:, ib] .= endnodeidxs[2, ib] + # teadjnodeidxs[ib, :] .= Int(endnodeidxs[2, ib]) + teadjnodeidxs[:, ib] .= Int(endnodeidxs[2, ib]) # teadjnodeidxs[2, ib] = -1 # set both dots and crosses to the same thing based on the single adjacent node. - tendotn[1, ib] = dot(tenormal[:, ib], normal[:, endpanelidxs[2, ib]]) + tendotn[1, ib] = dot(tenormal[:, ib], normal[:, Int(endpanelidxs[2, ib])]) tendotn[2, ib] = 0.0 # unnecessary since this value isn't used because the gamma value is prescribed to zero anyway, but just in case initialization changes - tencrossn[:, ib] .= cross2mag(normal[:, endpanelidxs[2, ib]], tenormal[:, ib]) + tencrossn[:, ib] .= cross2mag( + normal[:, Int(endpanelidxs[2, ib])], tenormal[:, ib] + ) end teinfluence_length[ib] = get_r(tenode[ib, 1, :], tenode[ib, 2, :])[2] end diff --git a/src/preprocess/geometry/rotor_geometry.jl b/src/preprocess/geometry/rotor_geometry.jl index bdb7aadd..d4074b0d 100644 --- a/src/preprocess/geometry/rotor_geometry.jl +++ b/src/preprocess/geometry/rotor_geometry.jl @@ -236,7 +236,7 @@ function generate_rotor_panels!( for irotor in 1:length(rotorzloc) @views xr[irotor] = [ fill(rotorzloc[irotor], nwake_sheets)' - wake_grid[2, rotor_indices_in_wake[irotor], 1:nwake_sheets]' + wake_grid[2, Int(rotor_indices_in_wake[irotor]), 1:nwake_sheets]' ] end diff --git a/src/preprocess/initialize_states.jl b/src/preprocess/initialize_states.jl new file mode 100644 index 00000000..deb482a0 --- /dev/null +++ b/src/preprocess/initialize_states.jl @@ -0,0 +1,457 @@ +""" +""" +function initialize_velocities( + operating_point, + blade_elements, + linsys, + ivr, + ivw, + body_totnodes, + wake_panel_sheet_be_map, +) + + ##### ----- Initialize ----- ##### + # - get type - # + #= + NOTE: use anything in the operating point, the wake on body AIC should cover any body and wake geometry changes, and the rotor-on-rotor velocity should cover any rotor changes. + =# + TF = promote_type( + eltype(operating_point.Omega), + eltype(operating_point.Vinf), + eltype(operating_point.rhoinf), + eltype(operating_point.muinf), + eltype(operating_point.asound), + eltype(linsys.A_bw), + eltype(ivr.v_rr), + eltype(blade_elements.twists), + eltype(blade_elements.chords), + eltype(blade_elements.Rtip), + eltype(blade_elements.Rhub), + ) + + # - Initialize velocity intermediates and outputs - # + # rename to only compute once + nbe = size(blade_elements.rotor_panel_centers) + nbe = size(blade_elements.rotor_panel_centers, 1) + + # outputs + vz_rotor = zeros(TF, nbe) + vtheta_rotor = zeros(TF, nbe) + Cm_wake = zeros(TF, size(ivw.v_ww, 1)) + + return initialize_velocities!( + vz_rotor, + vtheta_rotor, + Cm_wake, + operating_point, + blade_elements, + linsys, + ivr, + ivw, + body_totnodes, + wake_panel_sheet_be_map, + ) +end + +function initialize_velocities!( + vz_rotor, + vtheta_rotor, + Cm_wake, + operating_point, + blade_elements, + linsys, + ivr, + ivw, + body_totnodes, + wake_panel_sheet_be_map, +) + + # zero outputs: + vz_rotor .= 0 + vtheta_rotor .= 0 + Cm_wake .= 0 + + # - get floating point type - # + TF = promote_type( + eltype(operating_point.Omega), + eltype(operating_point.Vinf), + eltype(operating_point.rhoinf), + eltype(operating_point.muinf), + eltype(operating_point.asound), + eltype(linsys.A_bw), + eltype(ivr.v_rr), + eltype(blade_elements.twists), + eltype(blade_elements.chords), + eltype(blade_elements.Rtip), + eltype(blade_elements.Rhub), + ) + + # - Initialize intermediates - # + # rename to only compute once + nbe, nrotor = size(blade_elements.rotor_panel_centers) + + # intermediate values + # TODO: put these in a precomp container cache eventually + sigr = zeros(TF, nbe + 1, nrotor) + Cm_wake_vec = zeros(TF, nbe + 1) + vzind = zeros(TF, nbe) + vrind = zeros(TF, nbe) + vthetaind = zeros(TF, nbe) + + # Solve Linear System for gamb + # TODO; consider having an option here where you can fill the rhs cache (which should be used here) based on the reference velocity to try and get a better starting point + # #probably set that up in the precompute parameters function as this would be the first place that rhs vector would be seen. + gamb = ImplicitAD.implicit_linear( + linsys.A_bb, copy(linsys.b_bf); lsolve=ldiv!, Af=linsys.A_bb_LU + ) + + # - Get body-induced velocities on rotors - # + vzb = zeros(TF, nbe, nrotor) + vrb = zeros(TF, nbe, nrotor) + for irotor in 1:length(operating_point.Omega) + berange = (nbe * (irotor - 1) + 1):(nbe * irotor) + vzb[:, irotor] = ivr.v_rb[berange, :, 1] * gamb[1:body_totnodes] + vrb[:, irotor] = ivr.v_rb[berange, :, 2] * gamb[1:body_totnodes] + end + + ##### ----- Loop through rotors ----- ##### + for irotor in 1:length(operating_point.Omega) + + #remove body influece for previous rotor and add it for this rotor + if irotor > 1 + vzind .-= vzb[:, irotor - 1] + vrind .-= vrb[:, irotor - 1] + end + vzind .+= vzb[:, irotor] + vrind .+= vrb[:, irotor] + + # - Setup and Run CCBlade - # + # define rotor, do not apply any corrections (including a tip correction) + rotor = c4b.Rotor( + blade_elements.Rhub[irotor], + blade_elements.Rtip[irotor], + blade_elements.B[irotor]; + tip=nothing, + ) + + # define rotor sections + sections = + c4b.Section.( + blade_elements.rotor_panel_centers[:, irotor], + blade_elements.chords[:, irotor], + blade_elements.twists[:, irotor], + blade_elements.inner_airfoil[:, irotor], + ) + + # define operating points using induced velocity from rotors ahead of this one + c4bop = [ + c4b.OperatingPoint( + operating_point.Vinf[] + vz, # axial velocity V is freestream, vz is induced by bodies and rotor(s) ahead + operating_point.Omega[irotor] * + blade_elements.rotor_panel_centers[ir, irotor] + vt, # tangential velocity + operating_point.rhoinf[], + 0.0, #pitch is zero + operating_point.muinf[], + operating_point.asound[], + ) for (ir, (vz, vt)) in enumerate(zip(vzind, vthetaind)) + ] + + # solve CCBlade problem for this rotor + out = c4b.solve.(Ref(rotor), sections, c4bop) + + ##### ----- Assign Initial vz_rotor, vtheta_rotor, and Cm_wake ----- ##### + + # - vz_rotor and V_theta rotor - # + # self influence + vz_rotor[:, irotor] .+= vzind .+ getfield.(out, :u) + vtheta_rotor[:, irotor] .+= vthetaind .+ getfield.(out, :v) + + # - Get Cm_wake - # + #= + NOTE: we are going to estimate this by taking the velocities on the rotors (though using far field z terms) and applying them constantly straight back to the next rotor or end of wake. + =# + # Get source strengths + #= + NOTE: we need the values at the nodes not centers, so average the values and use the end values on the end points + =# + sigr[1] = @. blade_elements.B[irotor] / (4.0 * pi) * + getfield.(out, :cd)[1] * + getfield.(out, :W)[1] * + blade_elements.chords[1, irotor] + @. sigr[2:(end - 1)] = + ( + blade_elements.B[irotor] / (4.0 * pi) * + getfield.(out, :cd)[2:end] * + getfield.(out, :W)[2:end] * + blade_elements.chords[2:end, irotor] + + blade_elements.B[irotor] / (4.0 * pi) * + getfield.(out, :cd)[1:(end - 1)] * + getfield.(out, :W)[1:(end - 1)] * + blade_elements.chords[1:(end - 1), irotor] + ) / 2.0 + sigr[end] = @. blade_elements.B[irotor] / (4.0 * pi) * + getfield.(out, :cd)[end] * + getfield.(out, :W)[end] * + blade_elements.chords[end, irotor] + + # add influence of rotor radial induced velocity from self and rotors ahead + for jrotor in 1:irotor + berange = (nbe * (irotor - 1) + 1):(nbe * irotor) + berangej = ((nbe + 1) * (jrotor - 1) + 1):((nbe + 1) * jrotor) + vrind .+= ivr.v_rr[berange, berangej, 2] * sigr[:, jrotor] + end + + # add in axial and tangential influence aft of current rotor + vzind .+= 2.0 * getfield.(out, :u) + vthetaind .-= 2.0 * getfield.(out, :v) + + # since wakes extend from source panel endpoints, we need to average velocities and use the ends for endpoints + Cm_wake_vec[1] = sqrt((operating_point.Vinf[1] + vzind[1])^2 + vrind[1]^2) + Cm_wake_vec[2:(end - 1)] = + ( + sqrt.( + (operating_point.Vinf[1] .+ vzind[2:end]) .^ 2 .+ vrind[2:end] .^ 2 + ) .+ + sqrt.( + (operating_point.Vinf[1] .+ vzind[1:(end - 1)]) .^ 2 .+ + vrind[1:(end - 1)] .^ 2 + ) + ) / 2.0 + Cm_wake_vec[end] = sqrt((operating_point.Vinf[1] + vzind[end])^2 + vrind[end]^2) + + # fill in the section of the wake aft of the current rotor and up to the next rotor (or end of wake) + for (wid, wmap) in enumerate(eachrow(wake_panel_sheet_be_map)) + if wmap[2] >= irotor && wmap[2] < irotor + 1 + Cm_wake[wid] = Cm_wake_vec[wmap[1]] + end + end + end # loop through rotors + + return vz_rotor, vtheta_rotor, Cm_wake +end + +function initialize_strengths!( + Gamr, + sigr, + gamw, + operating_point, + blade_elements, + linsys, + ivr, + ivw, + wakeK, + body_totnodes, + wake_nodemap, + wake_endnodeidxs, + wake_panel_sheet_be_map, + wake_node_sheet_be_map, + wake_node_ids_along_casing_wake_interface, + wake_node_ids_along_centerbody_wake_interface, +) + + # zero outputs: + Gamr .= 0 + sigr .= 0 + gamw .= 0 + + # - get floating point type - # + TF = promote_type( + eltype(Gamr), + eltype(sigr), + eltype(gamw), + eltype(operating_point.Omega), + eltype(operating_point.Vinf), + eltype(operating_point.rhoinf), + eltype(operating_point.muinf), + eltype(operating_point.asound), + eltype(linsys.A_bw), + eltype(ivr.v_rr), + eltype(blade_elements.twists), + eltype(blade_elements.chords), + eltype(blade_elements.Rtip), + eltype(blade_elements.Rhub), + ) + + # - Initialize intermediates - # + # rename to only compute once + nbe, nrotor = size(blade_elements.rotor_panel_centers) + + # intermediate values + # TODO: put these in a precomp container cache eventually + vz_rotor = zeros(TF, nbe, nrotor) + vtheta_rotor = zeros(TF, nbe, nrotor) + Cm_wake_vec = zeros(TF, nbe + 1) + Cm_wake = zeros(TF, size(wake_panel_sheet_be_map, 1)) .= 0 + vthetaind = zeros(TF, nbe) + vzind = zeros(TF, nbe) + vrind = zeros(TF, nbe) + + # Solve Linear System for gamb + # 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) + + # - Get body-induced velocities on rotors - # + vzb = zeros(TF, nbe, nrotor) + vrb = zeros(TF, nbe, nrotor) + for irotor in 1:length(operating_point.Omega) + berange = (nbe * (irotor - 1) + 1):(nbe * irotor) + vzb[:, irotor] = ivr.v_rb[berange, :, 1] * gamb[1:body_totnodes] + vrb[:, irotor] = ivr.v_rb[berange, :, 2] * gamb[1:body_totnodes] + end + + ##### ----- Loop through rotors ----- ##### + for irotor in 1:nrotor + + #remove body influece for previous rotor and add it for this rotor + if irotor > 1 + vzind .-= vzb[:, irotor - 1] + vrind .-= vrb[:, irotor - 1] + end + vzind .+= vzb[:, irotor] + vrind .+= vrb[:, irotor] + + # - Setup and Run CCBlade - # + # define rotor, do not apply any corrections (including a tip correction) + rotor = c4b.Rotor( + blade_elements.Rhub[irotor], + blade_elements.Rtip[irotor], + blade_elements.B[irotor]; + tip=nothing, + ) + + # define rotor sections + sections = + c4b.Section.( + blade_elements.rotor_panel_centers[:, irotor], + blade_elements.chords[:, irotor], + blade_elements.twists[:, irotor], + blade_elements.inner_airfoil[:, irotor], + ) + + # define operating points using induced velocity from rotors ahead of this one + c4bop = [ + c4b.OperatingPoint( + operating_point.Vinf[] + vz, # axial velocity V is freestream, vz is induced by bodies and rotor(s) ahead + operating_point.Omega[irotor] * + blade_elements.rotor_panel_centers[ir, irotor] .- vt, # tangential velocity + operating_point.rhoinf[], + 0.0, #pitch is zero + operating_point.muinf[], + operating_point.asound[], + ) for (ir, (vz, vt)) in enumerate(zip(vzind, vthetaind)) + ] + + # solve CCBlade problem for this rotor + out = c4b.solve.(Ref(rotor), sections, c4bop) + + ##### ----- Assign Initial Gamr, sigr, and gamw ----- ##### + # - Get Gamr - # + Gamr[:, irotor] .= + 0.5 .* getfield.(out, :cl) .* getfield.(out, :W) .* + blade_elements.chords[:, irotor] + + # # - vz_rotor and V_theta rotor - # + # # self influence + # vz_rotor[:, irotor] .+= vzind .+ getfield.(out,:u + # vtheta_rotor[:, irotor] .+= vthetaind .+ getfield.(out,:v + + # - Get Cm_wake - # + #= + NOTE: we are going to estimate this by taking the velocities on the rotors (though using far field z terms) and applying them constantly straight back to the next rotor or end of wake. + =# + # Get source strengths + #= + NOTE: we need the values at the nodes not centers, so average the values and use the end values on the end points + =# + sigr[1, irotor] = @. blade_elements.B[irotor] / (4.0 * pi) * + getfield.(out, :cd)[1] * + getfield.(out, :W)[1] * + blade_elements.chords[1, irotor] + @. sigr[2:(end - 1), irotor] = + ( + blade_elements.B[irotor] / (4.0 * pi) * + getfield.(out, :cd)[2:end] * + getfield.(out, :W)[2:end] * + blade_elements.chords[2:end, irotor] + + blade_elements.B[irotor] / (4.0 * pi) * + getfield.(out, :cd)[1:(end - 1)] * + getfield.(out, :W)[1:(end - 1)] * + blade_elements.chords[1:(end - 1), irotor] + ) / 2.0 + sigr[end, irotor] = @. blade_elements.B[irotor] / (4.0 * pi) * + getfield.(out, :cd)[end] * + getfield.(out, :W)[end] * + blade_elements.chords[end, irotor] + + # add influence of rotor radial induced velocity from self and rotors ahead + for jrotor in 1:irotor + berange = (nbe * (irotor - 1) + 1):(nbe * irotor) + berangej = ((nbe + 1) * (jrotor - 1) + 1):((nbe + 1) * jrotor) + vrind .+= ivr.v_rr[berange, berangej, 2] * sigr[:, jrotor] + end + + # add in axial and tangential influence aft of current rotor + vzind .+= 2.0 * getfield.(out, :u) + vthetaind .-= 2.0 * getfield.(out, :v) + + # since wakes extend from source panel endpoints, we need to average velocities and use the ends for endpoints + Cm_wake_vec[1] = sqrt((operating_point.Vinf[1] + vzind[1])^2 + vrind[1]^2) + Cm_wake_vec[2:(end - 1)] = + ( + sqrt.( + (operating_point.Vinf[1] .+ vzind[2:end]) .^ 2 .+ vrind[2:end] .^ 2 + ) .+ + sqrt.( + (operating_point.Vinf[1] .+ vzind[1:(end - 1)]) .^ 2 .+ + vrind[1:(end - 1)] .^ 2 + ) + ) / 2.0 + Cm_wake_vec[end] = sqrt((operating_point.Vinf[1] + vzind[end])^2 + vrind[end]^2) + + # fill in the section of the wake aft of the current rotor and up to the next rotor (or end of wake) + for (wid, wmap) in enumerate(eachrow(wake_panel_sheet_be_map)) + if wmap[2] >= irotor && wmap[2] < irotor + 1 + Cm_wake[wid] = Cm_wake_vec[wmap[1]] + end + end + end # loop through rotors + + # - initialize wake strengths - # + # TODO: these should be in solve_containers, but need to figure out how to organize that as an input in this case + Gamma_tilde = zeros(TF, nbe, nrotor) + H_tilde = zeros(TF, nbe, nrotor) + deltaGamma2 = zeros(TF, nbe + 1, nrotor) + deltaH = zeros(TF, nbe + 1, nrotor) + Cm_avg = zeros(TF, size(gamw)) .= 0 + + average_wake_velocities!(Cm_avg, Cm_wake, wake_nodemap, wake_endnodeidxs) + + # - Calculate Wake Panel Strengths - # + # in-place solve for gamw, + calculate_wake_vortex_strengths!( + gamw, + Gamma_tilde, + H_tilde, + deltaGamma2, + deltaH, + Gamr, + Cm_avg, + blade_elements.B, + operating_point.Omega, + wakeK, + wake_node_sheet_be_map, + wake_node_ids_along_casing_wake_interface, + wake_node_ids_along_centerbody_wake_interface; + ) + + # Gamr struggles to converge if it's not initially positive... + for g in eachindex(Gamr) + if Gamr[g] < 0 + Gamr[g] = 0.05 + end + end + + return Gamr, sigr, gamw +end diff --git a/src/preprocess/preprocess.jl b/src/preprocess/preprocess.jl index 260e91bd..5d09ef4c 100644 --- a/src/preprocess/preprocess.jl +++ b/src/preprocess/preprocess.jl @@ -77,10 +77,6 @@ function reinterpolate_geometry!( paneling_constants; autoshiftduct=true, grid_solver_options=GridSolverOptions(), - # atol=1e-14, - # iteration_limit=100, - # relaxation_iteration_limit=3, - # relaxation_atol=1e-14, finterp=FLOWMath.akima, verbose=false, silence_warnings=true, @@ -156,10 +152,6 @@ function reinterpolate_geometry!( tip_gap[1], zwake; grid_solver_options=grid_solver_options, - # atol=atol, - # iteration_limit=iteration_limit, - # relaxation_iteration_limit=relaxation_iteration_limit, - # relaxation_atol=relaxation_atol, verbose=verbose, silence_warnings=silence_warnings, ) @@ -204,12 +196,12 @@ end """ function generate_all_panels!( panels, - idmaps, + wake_grid, rp_duct_coordinates, rp_centerbody_coordinates, - rotorstator_parameters, - paneling_constants, - wake_grid; + rotor_indices_in_wake, + rotorzloc, + nwake_sheets; itcpshift=0.05, axistol=1e-15, tegaptol=1e1 * eps(), @@ -218,8 +210,6 @@ function generate_all_panels!( # - Extract Tuples - # (; body_vortex_panels, rotor_source_panels, wake_vortex_panels) = panels - (; npanels, ncenterbody_inlet, nduct_inlet, wake_length, nwake_sheets, dte_minus_cbte) = - paneling_constants ##### ----- Fill Panel Objects ----- ##### # - Body Panels - # @@ -283,6 +273,11 @@ end """ """ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_options) + # - Reset Tuples - # + reset_containers!(ivr) + reset_containers!(ivw) + reset_containers!(ivb) + # - Extract Tuples - # # Extract induced velocities on rotor (; v_rr, v_rw, v_rb) = ivr @@ -307,7 +302,8 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_o body_vortex_panels.node, body_vortex_panels.nodemap, body_vortex_panels.influence_length, - ones(TF, 2, body_vortex_panels.totpanel), + #TODO: set things up differently so that you don't have to provide this allocated vector + ones(TF, 2, Int(body_vortex_panels.totpanel[])), integration_options, ) @@ -331,7 +327,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_o rotor_source_panels.node, rotor_source_panels.nodemap, rotor_source_panels.influence_length, - ones(TF, 2, rotor_source_panels.totnode), + ones(TF, 2, Int(rotor_source_panels.totnode[])), integration_options, ) @@ -343,7 +339,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_o wake_vortex_panels.node, wake_vortex_panels.nodemap, wake_vortex_panels.influence_length, - ones(TF, 2, wake_vortex_panels.totpanel), + ones(TF, 2, Int(wake_vortex_panels.totpanel[])), integration_options, ) @@ -368,7 +364,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_o rotor_source_panels.node, rotor_source_panels.nodemap, rotor_source_panels.influence_length, - ones(TF, 2, rotor_source_panels.totnode), + ones(TF, 2, Int(rotor_source_panels.totnode[])), integration_options, ) @@ -380,7 +376,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_o body_vortex_panels.node, body_vortex_panels.nodemap, body_vortex_panels.influence_length, - ones(TF, 2, body_vortex_panels.totpanel), + ones(TF, 2, Int(body_vortex_panels.totpanel[])), integration_options, ) @@ -404,7 +400,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_o wake_vortex_panels.node, wake_vortex_panels.nodemap, wake_vortex_panels.influence_length, - ones(TF, 2, wake_vortex_panels.totpanel), + ones(TF, 2, Int(wake_vortex_panels.totpanel[])), integration_options, ) @@ -430,7 +426,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_o body_vortex_panels.node, body_vortex_panels.nodemap, body_vortex_panels.influence_length, - ones(TF, 2, body_vortex_panels.totpanel), + ones(TF, 2, Int(body_vortex_panels.totpanel[])), integration_options, ) @@ -454,7 +450,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_o rotor_source_panels.node, rotor_source_panels.nodemap, rotor_source_panels.influence_length, - ones(TF, 2, rotor_source_panels.totpanel), + ones(TF, 2, Int(rotor_source_panels.totpanel[])), integration_options, ) @@ -466,7 +462,7 @@ function calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_o wake_vortex_panels.node, wake_vortex_panels.nodemap, wake_vortex_panels.influence_length, - ones(TF, 2, wake_vortex_panels.totpanel), + ones(TF, 2, Int(wake_vortex_panels.totpanel[])), integration_options, ) @@ -923,76 +919,6 @@ function interpolate_blade_elements!(blade_element_cache, rsp, rotor_panel_cente return (; outer_airfoil, inner_airfoil) end -# function interpolate_blade_elements!(blade_elements, rsp, Rtips, Rhubs, rotor_panel_center) -# # - Extract Blade Elements - # -# (; -# B, -# Rhub, -# Rtip, -# rotor_panel_centers, -# chords, -# twists, -# stagger, -# solidity, -# inner_fraction, -# fliplift, -# ) = blade_elements - -# Rtip .= Rtips -# Rhub .= Rhubs -# B .= rsp.B -# fliplift .= rsp.fliplift -# rotor_panel_centers .= rotor_panel_center - -# for irotor in 1:length(B) - -# # dimensionalize the blade element radial positions -# rblade = FLOWMath.linear([0.0; 1.0], [0.0; Rtip[irotor]], rsp.r) - -# # update chord lengths -# chords[:, irotor] .= FLOWMath.akima(rblade, chords, rotor_panel_centers) - -# # update twists -# twists[:, irotor] .= FLOWMath.akima(rblade, twists, rotor_panel_centers) - -# # update stagger -# stagger[:, irotor] .= get_stagger(twists) - -# # update solidity -# solidity[:, irotor] .= get_local_solidity(B, chords, rotor_panel_centers) - -# # get bounding airfoil polars -# outer_airfoil[:, irotor] .= similar(airfoils, length(rotor_panel_centers)) -# inner_airfoil[:, irotor] .= similar(airfoils, length(rotor_panel_centers)) -# inner_fraction[:, irotor] .= similar(airfoils, TF, length(rotor_panel_centers)) - -# for ir in 1:(length(rotor_panel_centers)) -# # outer airfoil -# io = min(length(rblade), searchsortedfirst(rblade, rotor_panel_centers[ir])) -# outer_airfoil[ir] = airfoils[io] - -# # inner airfoil -# ii = max(1, io - 1) -# inner_airfoil[ir] = airfoils[ii] - -# # fraction of inner airfoil's polars to use -# if rblade[io] == rblade[ii] -# inner_fraction[ir] = 1.0 -# else -# inner_fraction[ir] = -# (rotor_panel_centers[ir] - rblade[ii]) / (rblade[io] - rblade[ii]) -# end - -# # Check incorrect extrapolation -# if inner_fraction[ir] > 1.0 -# inner_fraction[ir] = 1.0 -# end -# end -# end - -# return (; inner_airfoil, outer_airfoil) -# end - """ wnm = wake_vortex_panels.nodemap wenids = wake_vortex_panels.endnodeidxs @@ -1165,7 +1091,7 @@ function set_index_maps( # - Map of wake node index to the wake sheet on which it resides and the last rotor ahead of the node - # nwsn = nwsp + 1 - wake_node_sheet_be_map = ones(Int, wenids[end], 2) + wake_node_sheet_be_map = ones(Int, Int(wenids[end]), 2) for i in 1:nwake_sheets wake_node_sheet_be_map[(1 + (i - 1) * nwsn):(i * nwsn), 1] .= i for (ir, r) in enumerate( @@ -1193,109 +1119,6 @@ function set_index_maps( ) end -""" -wnm = wake_vortex_panels.nodemap -venids = wake_vortex_panels.endnodeidxs -nbn = body_vortex_panels.totnode -wpsbm = wake_panel_sheet_be_map -e_map - wake_panel_sheet_be_map = ones(Int, wake_vortex_panels.totnode, 2) - num_wake_z_nodes = length(zwake) - for i in 1:(rotorstator_parameters[1].nwake_sheets) - wake_panel_sheet_be_map[(1 + (i - 1) * num_wake_z_nodes):(i * num_wake_z_nodes), 1] .= i - end - for (i, wn) in enumerate(eachcol(wake_vortex_panels.node)) - # TODO: DFDC geometry doesn't line up wake and rotor perfectly, so need a more robust option. - # wake_panel_sheet_be_map[i, 2] = findlast(x -> x <= wn[1], rotorstator_parameters.rotorzloc) - # TODO: current tee_mapsts are passing, but look here if things break in the future. - wake_panel_sheet_be_map[i, 2] = findmin(x -> abs(x - wn[1]), rotorstator_parameters.rotorzloc)[2] - end -""" -function set_index_maps!( - idmaps, - npanels, - nwake_sheets, - ncenterbody_inlet, - dte_minus_cbte, - wnm, - wenids, - nbn, - wpsbm, -) - # - Extract Index Maps - # - (; - wake_nodemap, - wake_endnodeidxs, - wake_panel_sheet_be_map, - wake_node_ids_along_casing_wake_interface, - wake_node_ids_along_centerbody_wake_interface, - id_of_first_casing_panel_aft_of_each_rotor, - id_of_first_centerbody_panel_aft_of_each_rotor, - rotor_indices_in_wake, - body_totnodes, - ) = idmaps - - wake_nodemap .= wnm - wake_endnodeidxs .= wenids - body_totnodes .= nbn - - if iszero(dte_minus_cbte) || dte_minus_cbte < 0 - cb_te_id = sum(npanels[1:(end - 1)]) + 1 - else - cb_te_id = sum(npanels[1:(end - 2)]) + 1 - end - wake_node_ids_along_centerbody_wake_interface .= collect(range(1, cb_te_id; step=1)) - - if iszero(dte_minus_cbte) || dte_minus_cbte > 0 - duct_te_id = sum(npanels[1:(end - 1)]) + 1 - else - duct_te_id = sum(npanels[1:(end - 2)]) + 1 - end - - ductteinwake = (sum(npanels) + 1) * nwake_sheets - (npanels[end] + 1) - - wake_node_ids_along_casing_wake_interface .= collect( - range(ductteinwake - duct_te_id + 1, ductteinwake; step=1) - ) - - for i in 1:nwake_sheets - wake_panel_sheet_be_map[(1 + (i - 1) * nwsn):(i * nwsn), 1] .= i - for (ir, r) in enumerate( - eachrow(@view(wake_panel_sheet_be_map[(1 + (i - 1) * nwsn):(i * nwsn), :])) - ) - r[2] = min(nrotor, searchsortedlast(rotor_indices_in_wake, ir)) - end - end - - if dte_minus_cbte < 0 - id_of_first_casing_panel_aft_of_each_rotor = cumsum([ - npanels[i] for i in (length(npanels) - 1):-1:1 - ])[2:end] - elseif dte_minus_cbte > 0 - id_of_first_casing_panel_aft_of_each_rotor = cumsum([ - npanels[i] for i in (length(npanels) - 2):-1:1 - ]) - else - id_of_first_casing_panel_aft_of_each_rotor = cumsum([ - npanels[i] for i in (length(npanels) - 1):-1:1 - ]) - end - - if iszero(dte_minus_cbte) - id_of_first_centerbody_panel_aft_of_each_rotor = [ - ncenterbody_inlet + 1 - ncenterbody_inlet .+ [npanels[i] for i in 1:(length(npanels) - 2)] - ] - else - id_of_first_centerbody_panel_aft_of_each_rotor = [ - ncenterbody_inlet + 1 - ncenterbody_inlet .+ [npanels[i] for i in 1:(length(npanels) - 3)] - ] - end - - return idmaps -end - """ """ function precompute_parameters( @@ -1487,7 +1310,9 @@ function precompute_parameters!( blade_element_cache, linsys, wakeK, - propulsor; + propulsor, + prepost_containers, + problem_dimensions; grid_solver_options=GridSolverOptions(), integration_options=IntegrationOptions(), autoshiftduct=true, @@ -1499,7 +1324,7 @@ function precompute_parameters!( verbose=false, ) - # - Extract propulsor - # + # - unpack propulsor - # (; duct_coordinates, # Matrix centerbody_coordinates, # Matrix @@ -1509,18 +1334,27 @@ function precompute_parameters!( reference_parameters, # NamedTuple of vectors (TODO) of numbers ) = propulsor - # - Get Problem Dimensions - # - problem_dimensions = get_problem_dimensions(paneling_constants) + # - Unpack preprocess containers - # + (; + wake_grid, + rp_duct_coordinates, + rp_centerbody_coordinates, + rotor_indices_in_wake, + ) = prepost_containers # - Reinterpolate Geometry and Generate Wake Grid - # - wake_grid, rp_duct_coordinates, rp_centerbody_coordinates, rotor_indices_in_wake = reinterpolate_geometry( - problem_dimensions, + # note: currently has 3526 allocations and takes 103.21 ms + reinterpolate_geometry!( + wake_grid, + rp_duct_coordinates, + rp_centerbody_coordinates, + rotor_indices_in_wake, duct_coordinates, centerbody_coordinates, rotorstator_parameters, paneling_constants; - grid_solver_options=grid_solver_options, autoshiftduct=autoshiftduct, + grid_solver_options=grid_solver_options, finterp=finterp, verbose=verbose, silence_warnings=silence_warnings, @@ -1550,6 +1384,7 @@ function precompute_parameters!( paneling_constants, operating_point, reference_parameters, + prepost_containers, problem_dimensions; grid_solver_options=grid_solver_options, integration_options=integration_options, @@ -1579,6 +1414,7 @@ function precompute_parameters!( paneling_constants, operating_point, reference_parameters, + prepost_containers, problem_dimensions=nothing; grid_solver_options=GridSolverOptions(), integration_options=IntegrationOptions(), @@ -1597,6 +1433,10 @@ function precompute_parameters!( reset_containers!(blade_element_cache) reset_containers!(linsys) reset_containers!(wakeK) + reset_containers!(prepost_containers; exception_keys=[:wake_grid, + :rp_duct_coordinates, + :rp_centerbody_coordinates, + :rotor_indices_in_wake]) # - Get Floating Point Type - # TF = promote_type( @@ -1615,14 +1455,25 @@ function precompute_parameters!( eltype(rotorstator_parameters.twists), ) + # - Unpack preprocess containers - # + (; + panels, + ivb, + AICn, + AICpcp, + vdnb, + vdnpcp, + ) = prepost_containers + # - Panel Everything - # - body_vortex_panels, rotor_source_panels, wake_vortex_panels = generate_all_panels( + generate_all_panels!( + panels, + wake_grid, rp_duct_coordinates, rp_centerbody_coordinates, - paneling_constants.nwake_sheets, rotor_indices_in_wake, rotorstator_parameters.rotorzloc, - wake_grid; + paneling_constants.nwake_sheets; itcpshift=itcpshift, axistol=axistol, tegaptol=tegaptol, @@ -1632,62 +1483,35 @@ function precompute_parameters!( # - Get problem dimensions if not already done - # if isnothing(problem_dimensions) problem_dimensions = get_problem_dimensions( - body_vortex_panels, rotor_source_panels, wake_vortex_panels + panels.body_vortex_panels, panels.rotor_source_panels, panels.wake_vortex_panels ) end # - Compute Influence Matrices - # - # TODO: put ivb in post-process cache eventually - ivb = (; - v_bb=zeros(TF, problem_dimensions.nbp, problem_dimensions.nbn, 2), - v_br=zeros( - TF, - problem_dimensions.nbp, - problem_dimensions.nrotor * problem_dimensions.nws, - 2, - ), - 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), - integration_options, - ) + calculate_unit_induced_velocities!(ivr, ivw, ivb, panels, integration_options) # - Set up Linear System - # - # TODO: put these containers in a precomp cache eventually - intermediate_containers = (; - AICn=zeros(TF, problem_dimensions.nbp, problem_dimensions.nbn), - AICpcp=zeros(TF, 2, problem_dimensions.nbn), - vdnb=zeros(TF, problem_dimensions.nbp), - vdnpcp=zeros(TF, 2), - ) - - # TODO: test this function A_bb_LU, lu_decomp_flag = initialize_linear_system!( linsys, ivb, - body_vortex_panels, - rotor_source_panels, - wake_vortex_panels, + panels.body_vortex_panels, + panels.rotor_source_panels, + panels.wake_vortex_panels, operating_point.Vinf[1], - intermediate_containers, - integration_options + (; AICn, AICpcp, vdnb, vdnpcp), + integration_options, ) # - Interpolate Blade Elements - # - # TODO: test this function airfoils = interpolate_blade_elements!( blade_element_cache, rotorstator_parameters, - rotor_source_panels.controlpoint[2, :], + panels.rotor_source_panels.controlpoint[2, :], problem_dimensions.nbe, ) # - Get geometry-based constants for wake node strength calculations - # - get_wake_k!(wakeK, wake_vortex_panels.node[2, :]) + get_wake_k!(wakeK, panels.wake_vortex_panels.node[2, :]) # - Save all the index mapping (bookkeeping) - # idmaps = set_index_maps( @@ -1695,8 +1519,8 @@ function precompute_parameters!( paneling_constants.ncenterbody_inlet, paneling_constants.nwake_sheets, paneling_constants.dte_minus_cbte, - wake_vortex_panels.nodemap, - wake_vortex_panels.endnodeidxs, + Int.(panels.wake_vortex_panels.nodemap), + Int.(panels.wake_vortex_panels.endnodeidxs), problem_dimensions.nwp, problem_dimensions.nwsp, problem_dimensions.nbn, @@ -1705,559 +1529,5 @@ function precompute_parameters!( problem_dimensions.nrotor, ) - return ivb, - A_bb_LU, - lu_decomp_flag, - airfoils, - idmaps, - (; body_vortex_panels, rotor_source_panels, wake_vortex_panels), - problem_dimensions -end - -# """ -# """ -# function TODOprecompute_parameters!( -# ivr, -# ivw, -# ivb, -# linsys, -# blade_elements, -# idmaps, -# panels, -# propulsor, -# precomp_containers; # contains wake_grid and repaneled duct and centerbody coordinates -# grid_solver_options=GridSolverOptions(), -# itcpshift=0.05, -# axistol=1e-15, -# tegaptol=1e1 * eps(), -# finterp=FLOWMath.akima, -# silence_warnings=true, -# verbose=false, -# ) - -# # - Extract propulsor - # -# (; -# duct_coordinates, # Matrix -# centerbody_coordinates, # Matrix -# rotorstator_parameters, # NamedTuple of vectors of numbers -# paneling_constants, # NamedTuple of scalars and vectors of scalars -# operating_point, # NamedTuple of vectors of numbers -# reference_parameters, # NamedTuple of vectors of numbers -# ) = propulsor - -# # - Extract Precomp Containers - # -# (; panels, wake_grid, rp_duct_coordinates, rp_centerbody_coordinates, AICn, AICpcp) = -# precomp_containers - -# # - Reinterpolate Geometry and Generate Wake Grid - # -# # TODO: test this function -# reinterpolate_geometry!( -# precomp_containers, -# duct_coordinates, -# centerbody_coordinates, -# rotorstator_parameters, -# paneling_constants, -# idmaps.rotor_indices_in_wake; -# grid_solver_options=grid_solver_options, -# finterp=finterp, -# silence_warnings=silence_warnings, -# ) - -# # - Panel Everything - # -# # TODO: test this function -# generate_all_panels!( -# panels, -# idmaps, -# rp_duct_coordinates, -# rp_centerbody_coordinates, -# rotorstator_parameters, -# paneling_constants, -# wake_grid; -# itcpshift=itcpshift, -# axistol=axistol, -# tegaptol=tegaptol, -# silence_warnings=silence_warnings, -# ) - -# # - Compute Influence Matrices - # -# # TODO: test this function -# 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, integration_options -# ) - -# # - Interpolate Blade Elements - # -# # TODO: test this function -# airfoils = interpolate_blade_elements!(blade_elements, panels.rotor_source_panels) - -# # - Save all the index mapping (bookkeeping) - # -# # TODO: test this function -# set_index_maps!(idmaps, paneling_constants) - -# # - Get geometry-based constants for wake node strength calculations - # -# # TODO: write in-place version of this function -# get_wake_k!(wakeK, wake_vortex_panels.node[2, :]) - -# return ivr, ivw, ivb, linsys, blade_elements, wakeK, idmaps, panels, lu_decomp_flag -# end - -""" -""" -function initialize_velocities( - operating_point, - blade_elements, - linsys, - ivr, - ivw, - body_totnodes, - wake_panel_sheet_be_map, -) - - ##### ----- Initialize ----- ##### - # - get type - # - #= - NOTE: use anything in the operating point, the wake on body AIC should cover any body and wake geometry changes, and the rotor-on-rotor velocity should cover any rotor changes. - =# - TF = promote_type( - eltype(operating_point.Omega), - eltype(operating_point.Vinf), - eltype(operating_point.rhoinf), - eltype(operating_point.muinf), - eltype(operating_point.asound), - eltype(linsys.A_bw), - eltype(ivr.v_rr), - eltype(blade_elements.twists), - eltype(blade_elements.chords), - eltype(blade_elements.Rtip), - eltype(blade_elements.Rhub), - ) - - # - Initialize velocity intermediates and outputs - # - # rename to only compute once - nbe = size(blade_elements.rotor_panel_centers) - nbe = size(blade_elements.rotor_panel_centers, 1) - - # outputs - vz_rotor = zeros(TF, nbe) - vtheta_rotor = zeros(TF, nbe) - Cm_wake = zeros(TF, size(ivw.v_ww, 1)) - - return initialize_velocities!( - vz_rotor, - vtheta_rotor, - Cm_wake, - operating_point, - blade_elements, - linsys, - ivr, - ivw, - body_totnodes, - wake_panel_sheet_be_map, - ) -end - -function initialize_velocities!( - vz_rotor, - vtheta_rotor, - Cm_wake, - operating_point, - blade_elements, - linsys, - ivr, - ivw, - body_totnodes, - wake_panel_sheet_be_map, -) - - # zero outputs: - vz_rotor .= 0 - vtheta_rotor .= 0 - Cm_wake .= 0 - - # - get floating point type - # - TF = promote_type( - eltype(operating_point.Omega), - eltype(operating_point.Vinf), - eltype(operating_point.rhoinf), - eltype(operating_point.muinf), - eltype(operating_point.asound), - eltype(linsys.A_bw), - eltype(ivr.v_rr), - eltype(blade_elements.twists), - eltype(blade_elements.chords), - eltype(blade_elements.Rtip), - eltype(blade_elements.Rhub), - ) - - # - Initialize intermediates - # - # rename to only compute once - nbe, nrotor = size(blade_elements.rotor_panel_centers) - - # intermediate values - # TODO: put these in a precomp container cache eventually - sigr = zeros(TF, nbe + 1, nrotor) - Cm_wake_vec = zeros(TF, nbe + 1) - vzind = zeros(TF, nbe) - vrind = zeros(TF, nbe) - vthetaind = zeros(TF, nbe) - - # Solve Linear System for gamb - # TODO; consider having an option here where you can fill the rhs cache (which should be used here) based on the reference velocity to try and get a better starting point - # #probably set that up in the precompute parameters function as this would be the first place that rhs vector would be seen. - gamb = ImplicitAD.implicit_linear( - linsys.A_bb, copy(linsys.b_bf); lsolve=ldiv!, Af=linsys.A_bb_LU - ) - - # - Get body-induced velocities on rotors - # - vzb = zeros(TF, nbe, nrotor) - vrb = zeros(TF, nbe, nrotor) - for irotor in 1:length(operating_point.Omega) - berange = (nbe * (irotor - 1) + 1):(nbe * irotor) - vzb[:, irotor] = ivr.v_rb[berange, :, 1] * gamb[1:body_totnodes] - vrb[:, irotor] = ivr.v_rb[berange, :, 2] * gamb[1:body_totnodes] - end - - ##### ----- Loop through rotors ----- ##### - for irotor in 1:length(operating_point.Omega) - - #remove body influece for previous rotor and add it for this rotor - if irotor > 1 - vzind .-= vzb[:, irotor - 1] - vrind .-= vrb[:, irotor - 1] - end - vzind .+= vzb[:, irotor] - vrind .+= vrb[:, irotor] - - # - Setup and Run CCBlade - # - # define rotor, do not apply any corrections (including a tip correction) - rotor = c4b.Rotor( - blade_elements.Rhub[irotor], - blade_elements.Rtip[irotor], - blade_elements.B[irotor]; - tip=nothing, - ) - - # define rotor sections - sections = - c4b.Section.( - blade_elements.rotor_panel_centers[:, irotor], - blade_elements.chords[:, irotor], - blade_elements.twists[:, irotor], - blade_elements.inner_airfoil[:, irotor], - ) - - # define operating points using induced velocity from rotors ahead of this one - c4bop = [ - c4b.OperatingPoint( - operating_point.Vinf[] + vz, # axial velocity V is freestream, vz is induced by bodies and rotor(s) ahead - operating_point.Omega[irotor] * - blade_elements.rotor_panel_centers[ir, irotor] + vt, # tangential velocity - operating_point.rhoinf[], - 0.0, #pitch is zero - operating_point.muinf[], - operating_point.asound[], - ) for (ir, (vz, vt)) in enumerate(zip(vzind, vthetaind)) - ] - - # solve CCBlade problem for this rotor - out = c4b.solve.(Ref(rotor), sections, c4bop) - - ##### ----- Assign Initial vz_rotor, vtheta_rotor, and Cm_wake ----- ##### - - # - vz_rotor and V_theta rotor - # - # self influence - vz_rotor[:, irotor] .+= vzind .+ getfield.(out, :u) - vtheta_rotor[:, irotor] .+= vthetaind .+ getfield.(out, :v) - - # - Get Cm_wake - # - #= - NOTE: we are going to estimate this by taking the velocities on the rotors (though using far field z terms) and applying them constantly straight back to the next rotor or end of wake. - =# - # Get source strengths - #= - NOTE: we need the values at the nodes not centers, so average the values and use the end values on the end points - =# - sigr[1] = @. blade_elements.B[irotor] / (4.0 * pi) * - getfield.(out, :cd)[1] * - getfield.(out, :W)[1] * - blade_elements.chords[1, irotor] - @. sigr[2:(end - 1)] = - ( - blade_elements.B[irotor] / (4.0 * pi) * - getfield.(out, :cd)[2:end] * - getfield.(out, :W)[2:end] * - blade_elements.chords[2:end, irotor] + - blade_elements.B[irotor] / (4.0 * pi) * - getfield.(out, :cd)[1:(end - 1)] * - getfield.(out, :W)[1:(end - 1)] * - blade_elements.chords[1:(end - 1), irotor] - ) / 2.0 - sigr[end] = @. blade_elements.B[irotor] / (4.0 * pi) * - getfield.(out, :cd)[end] * - getfield.(out, :W)[end] * - blade_elements.chords[end, irotor] - - # add influence of rotor radial induced velocity from self and rotors ahead - for jrotor in 1:irotor - berange = (nbe * (irotor - 1) + 1):(nbe * irotor) - berangej = ((nbe + 1) * (jrotor - 1) + 1):((nbe + 1) * jrotor) - vrind .+= ivr.v_rr[berange, berangej, 2] * sigr[:, jrotor] - end - - # add in axial and tangential influence aft of current rotor - vzind .+= 2.0 * getfield.(out, :u) - vthetaind .-= 2.0 * getfield.(out, :v) - - # since wakes extend from source panel endpoints, we need to average velocities and use the ends for endpoints - Cm_wake_vec[1] = sqrt((operating_point.Vinf[1] + vzind[1])^2 + vrind[1]^2) - Cm_wake_vec[2:(end - 1)] = - ( - sqrt.( - (operating_point.Vinf[1] .+ vzind[2:end]) .^ 2 .+ vrind[2:end] .^ 2 - ) .+ - sqrt.( - (operating_point.Vinf[1] .+ vzind[1:(end - 1)]) .^ 2 .+ - vrind[1:(end - 1)] .^ 2 - ) - ) / 2.0 - Cm_wake_vec[end] = sqrt((operating_point.Vinf[1] + vzind[end])^2 + vrind[end]^2) - - # fill in the section of the wake aft of the current rotor and up to the next rotor (or end of wake) - for (wid, wmap) in enumerate(eachrow(wake_panel_sheet_be_map)) - if wmap[2] >= irotor && wmap[2] < irotor + 1 - Cm_wake[wid] = Cm_wake_vec[wmap[1]] - end - end - end # loop through rotors - - return vz_rotor, vtheta_rotor, Cm_wake -end - -function initialize_strengths!( - Gamr, - sigr, - gamw, - operating_point, - blade_elements, - linsys, - ivr, - ivw, - wakeK, - body_totnodes, - wake_nodemap, - wake_endnodeidxs, - wake_panel_sheet_be_map, - wake_node_sheet_be_map, - wake_node_ids_along_casing_wake_interface, - wake_node_ids_along_centerbody_wake_interface, -) - - # zero outputs: - Gamr .= 0 - sigr .= 0 - gamw .= 0 - - # - get floating point type - # - TF = promote_type( - eltype(Gamr), - eltype(sigr), - eltype(gamw), - eltype(operating_point.Omega), - eltype(operating_point.Vinf), - eltype(operating_point.rhoinf), - eltype(operating_point.muinf), - eltype(operating_point.asound), - eltype(linsys.A_bw), - eltype(ivr.v_rr), - eltype(blade_elements.twists), - eltype(blade_elements.chords), - eltype(blade_elements.Rtip), - eltype(blade_elements.Rhub), - ) - - # - Initialize intermediates - # - # rename to only compute once - nbe, nrotor = size(blade_elements.rotor_panel_centers) - - # intermediate values - # TODO: put these in a precomp container cache eventually - vz_rotor = zeros(TF, nbe, nrotor) - vtheta_rotor = zeros(TF, nbe, nrotor) - Cm_wake_vec = zeros(TF, nbe + 1) - Cm_wake = zeros(TF, size(wake_panel_sheet_be_map, 1)) .= 0 - vthetaind = zeros(TF, nbe) - vzind = zeros(TF, nbe) - vrind = zeros(TF, nbe) - - # Solve Linear System for gamb - # 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) - - # - Get body-induced velocities on rotors - # - vzb = zeros(TF, nbe, nrotor) - vrb = zeros(TF, nbe, nrotor) - for irotor in 1:length(operating_point.Omega) - berange = (nbe * (irotor - 1) + 1):(nbe * irotor) - vzb[:, irotor] = ivr.v_rb[berange, :, 1] * gamb[1:body_totnodes] - vrb[:, irotor] = ivr.v_rb[berange, :, 2] * gamb[1:body_totnodes] - end - - ##### ----- Loop through rotors ----- ##### - for irotor in 1:nrotor - - #remove body influece for previous rotor and add it for this rotor - if irotor > 1 - vzind .-= vzb[:, irotor - 1] - vrind .-= vrb[:, irotor - 1] - end - vzind .+= vzb[:, irotor] - vrind .+= vrb[:, irotor] - - # - Setup and Run CCBlade - # - # define rotor, do not apply any corrections (including a tip correction) - rotor = c4b.Rotor( - blade_elements.Rhub[irotor], - blade_elements.Rtip[irotor], - blade_elements.B[irotor]; - tip=nothing, - ) - - # define rotor sections - sections = - c4b.Section.( - blade_elements.rotor_panel_centers[:, irotor], - blade_elements.chords[:, irotor], - blade_elements.twists[:, irotor], - blade_elements.inner_airfoil[:, irotor], - ) - - # define operating points using induced velocity from rotors ahead of this one - c4bop = [ - c4b.OperatingPoint( - operating_point.Vinf[] + vz, # axial velocity V is freestream, vz is induced by bodies and rotor(s) ahead - operating_point.Omega[irotor] * - blade_elements.rotor_panel_centers[ir, irotor] .- vt, # tangential velocity - operating_point.rhoinf[], - 0.0, #pitch is zero - operating_point.muinf[], - operating_point.asound[], - ) for (ir, (vz, vt)) in enumerate(zip(vzind, vthetaind)) - ] - - # solve CCBlade problem for this rotor - out = c4b.solve.(Ref(rotor), sections, c4bop) - - ##### ----- Assign Initial Gamr, sigr, and gamw ----- ##### - # - Get Gamr - # - Gamr[:, irotor] .= - 0.5 .* getfield.(out, :cl) .* getfield.(out, :W) .* - blade_elements.chords[:, irotor] - - # # - vz_rotor and V_theta rotor - # - # # self influence - # vz_rotor[:, irotor] .+= vzind .+ getfield.(out,:u - # vtheta_rotor[:, irotor] .+= vthetaind .+ getfield.(out,:v - - # - Get Cm_wake - # - #= - NOTE: we are going to estimate this by taking the velocities on the rotors (though using far field z terms) and applying them constantly straight back to the next rotor or end of wake. - =# - # Get source strengths - #= - NOTE: we need the values at the nodes not centers, so average the values and use the end values on the end points - =# - sigr[1, irotor] = @. blade_elements.B[irotor] / (4.0 * pi) * - getfield.(out, :cd)[1] * - getfield.(out, :W)[1] * - blade_elements.chords[1, irotor] - @. sigr[2:(end - 1), irotor] = - ( - blade_elements.B[irotor] / (4.0 * pi) * - getfield.(out, :cd)[2:end] * - getfield.(out, :W)[2:end] * - blade_elements.chords[2:end, irotor] + - blade_elements.B[irotor] / (4.0 * pi) * - getfield.(out, :cd)[1:(end - 1)] * - getfield.(out, :W)[1:(end - 1)] * - blade_elements.chords[1:(end - 1), irotor] - ) / 2.0 - sigr[end, irotor] = @. blade_elements.B[irotor] / (4.0 * pi) * - getfield.(out, :cd)[end] * - getfield.(out, :W)[end] * - blade_elements.chords[end, irotor] - - # add influence of rotor radial induced velocity from self and rotors ahead - for jrotor in 1:irotor - berange = (nbe * (irotor - 1) + 1):(nbe * irotor) - berangej = ((nbe + 1) * (jrotor - 1) + 1):((nbe + 1) * jrotor) - vrind .+= ivr.v_rr[berange, berangej, 2] * sigr[:, jrotor] - end - - # add in axial and tangential influence aft of current rotor - vzind .+= 2.0 * getfield.(out, :u) - vthetaind .-= 2.0 * getfield.(out, :v) - - # since wakes extend from source panel endpoints, we need to average velocities and use the ends for endpoints - Cm_wake_vec[1] = sqrt((operating_point.Vinf[1] + vzind[1])^2 + vrind[1]^2) - Cm_wake_vec[2:(end - 1)] = - ( - sqrt.( - (operating_point.Vinf[1] .+ vzind[2:end]) .^ 2 .+ vrind[2:end] .^ 2 - ) .+ - sqrt.( - (operating_point.Vinf[1] .+ vzind[1:(end - 1)]) .^ 2 .+ - vrind[1:(end - 1)] .^ 2 - ) - ) / 2.0 - Cm_wake_vec[end] = sqrt((operating_point.Vinf[1] + vzind[end])^2 + vrind[end]^2) - - # fill in the section of the wake aft of the current rotor and up to the next rotor (or end of wake) - for (wid, wmap) in enumerate(eachrow(wake_panel_sheet_be_map)) - if wmap[2] >= irotor && wmap[2] < irotor + 1 - Cm_wake[wid] = Cm_wake_vec[wmap[1]] - end - end - end # loop through rotors - - # - initialize wake strengths - # - # TODO: these should be in solve_containers, but need to figure out how to organize that as an input in this case - Gamma_tilde = zeros(TF, nbe, nrotor) - H_tilde = zeros(TF, nbe, nrotor) - deltaGamma2 = zeros(TF, nbe + 1, nrotor) - deltaH = zeros(TF, nbe + 1, nrotor) - Cm_avg = zeros(TF, size(gamw)) .= 0 - - average_wake_velocities!(Cm_avg, Cm_wake, wake_nodemap, wake_endnodeidxs) - - # - Calculate Wake Panel Strengths - # - # in-place solve for gamw, - calculate_wake_vortex_strengths!( - gamw, - Gamma_tilde, - H_tilde, - deltaGamma2, - deltaH, - Gamr, - Cm_avg, - blade_elements.B, - operating_point.Omega, - wakeK, - wake_node_sheet_be_map, - wake_node_ids_along_casing_wake_interface, - wake_node_ids_along_centerbody_wake_interface; - ) - - # Gamr struggles to converge if it's not initially positive... - for g in eachindex(Gamr) - if Gamr[g] < 0 - Gamr[g] = 0.05 - end - end - - return Gamr, sigr, gamw + return ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps, problem_dimensions end diff --git a/src/preprocess/velocities/body_aic.jl b/src/preprocess/velocities/body_aic.jl index 53f90079..090a2b2b 100644 --- a/src/preprocess/velocities/body_aic.jl +++ b/src/preprocess/velocities/body_aic.jl @@ -86,8 +86,8 @@ function vortex_aic_boundary_on_boundary!( for (j, (nmap, lj)) in enumerate(zip(eachcol(nodemap), influence_length)) # Loop through control points being influenced for (i, (cpi, nhat)) in enumerate(zip(eachcol(controlpoint), eachcol(normal))) - n1 = view(node, :, nmap[1]) - n2 = view(node, :, nmap[2]) + n1 = view(node, :, Int(nmap[1])) + n2 = view(node, :, Int(nmap[2])) # get unit induced velocity from the panel onto the control point if i != j @@ -118,7 +118,7 @@ function vortex_aic_boundary_on_boundary!( for k in 1:2 # fill the Matrix - AICn[i, nmap[k]] += dot(vel[k, :], nhat) + AICn[i, Int(nmap[k])] += dot(vel[k, :], nhat) end #for k end #for i end #for j @@ -209,8 +209,8 @@ function vortex_aic_boundary_on_field!( for (i, (cpi, nhat)) in enumerate(zip(eachcol(controlpoint), eachcol(normal))) # loop through panels doing the influencing for (j, (nmap, lj)) in enumerate(zip(eachcol(nodemap), influence_length)) - n1 = view(node, :, nmap[1]) - n2 = view(node, :, nmap[2]) + n1 = view(node, :, Int(nmap[1])) + n2 = view(node, :, Int(nmap[2])) # check of self-induced: if isapprox(cpi, 0.5 * (n1 .+ n2)) @@ -249,7 +249,7 @@ function vortex_aic_boundary_on_field!( for k in 1:2 # fill the Matrix - AICn[i, nmap[k]] += dot(vel[k, :], nhat) + AICn[i, Int(nmap[k])] += dot(vel[k, :], nhat) end #for k end #for j end #for i @@ -343,11 +343,11 @@ function add_te_gap_aic!( for k in 1:2 # fill the Matrix - AICn[i, nmap[k]] += dot(ndn[k] * vvel[k, :] + ncn[k] * svel[k, :], nhat) + AICn[i, Int(nmap[k])] += dot(ndn[k] * vvel[k, :] + ncn[k] * svel[k, :], nhat) #TODO: is this a bug? seems like the above line should be in an else statement if wake # wake "TE Panels" only have the vortex influence - AICn[i, nmap[k]] += dot(ndn[k] * vvel[k, :], nhat) + AICn[i, Int(nmap[k])] += dot(ndn[k] * vvel[k, :], nhat) end end #for k end #for j @@ -439,8 +439,8 @@ function source_aic!( for (i, (cpi, nhat)) in enumerate(zip(eachcol(controlpoint), eachcol(normal))) # loop through panels doing the influencing for (j, (nmap, lj)) in enumerate(zip(eachcol(nodemap), influence_length)) - n1 = view(node, :, nmap[1]) - n2 = view(node, :, nmap[2]) + n1 = view(node, :, Int(nmap[1])) + n2 = view(node, :, Int(nmap[2])) # get unit induced velocity from the panel onto the control point # vel .= nominal_vortex_panel_integration(n1, n2, lj, cpi) @@ -452,7 +452,7 @@ function source_aic!( for k in 1:2 # fill the Matrix - AICn[i, nmap[k]] += dot(vel[k, :], nhat) + AICn[i, Int(nmap[k])] += dot(vel[k, :], nhat) end #for k end #for j end #for i @@ -541,30 +541,30 @@ function assemble_lhs_matrix!( # - Place standard influence coefficients - # # 1:totpanel, 1:totnode are standard AIC values - LHS[1:totpanel, 1:totnode] .= AICn + LHS[1:Int(totpanel[]), 1:Int(totnode[])] .= AICn # - Place Dummy influences on right-most columns - # # Duct Dummy's - LHS[1:npanel[1], end - 1] .= dummyval + LHS[1:Int(npanel[1]), end - 1] .= dummyval # Hub Dummy's - LHS[(npanel[1] + 1):totpanel, end] .= dummyval + LHS[(Int(npanel[1]) + 1):Int(totpanel[]), end] .= dummyval # - Place internal duct pseudo control point row - # - LHS[totpanel + 1, 1:totnode] .= @view(AICpcp[1, :]) + LHS[Int(totpanel[]) + 1, 1:Int(totnode[])] .= @view(AICpcp[1, :]) # - Place Kutta condition Row - # - LHS[(totpanel + 2), 1] = LHS[totpanel + 2, nnode[1]] = 1.0 + LHS[(Int(totpanel[]) + 2), 1] = LHS[Int(totpanel[]) + 2, Int(nnode[1])] = 1.0 # - Place hub LE prescribed panel row - # - LHS[totpanel + 3, prescribednodeidxs[1]] = 1.0 + LHS[Int(totpanel[]) + 3, Int(prescribednodeidxs[1])] = 1.0 # - Place hub TE prescribed panel row OR hub internal pseudo control point row - # - if length(prescribednodeidxs) > 1 - # prescribed TE node - LHS[totpanel + 4, prescribednodeidxs[2]] = 1.0 - else + if iszero(Int(prescribednodeidxs[2])) # internal pseudo control point - LHS[totpanel + 4, 1:totnode] .= @view(AICpcp[2, :]) + LHS[Int(totpanel[]) + 4, 1:Int(totnode[])] .= @view(AICpcp[2, :]) + else + # prescribed TE node + LHS[Int(totpanel[]) + 4, Int(prescribednodeidxs[2])] = 1.0 end return LHS @@ -637,14 +637,14 @@ function assemble_rhs_matrix!( # - Place standard influence coefficients - # # 1:totpanel, 1:totnode are standard AIC values - RHS[1:totpanel] .= vdnb + RHS[1:Int(totpanel[])] .= vdnb # - Place Duct pseudo control point freestream influence - # - RHS[totnode - 1] = vdnpcp[1] + RHS[Int(totnode[]) - 1] = vdnpcp[1] # - Place hub pseudo control point freestream influence - # - if length(prescribednodeidxs) == 1 - RHS[totnode + 2] = vdnpcp[2] + if iszero(Int(prescribednodeidxs[2])) + RHS[Int(totnode[]) + 2] = vdnpcp[2] end return RHS @@ -664,7 +664,7 @@ function assemble_rhs_matrix!( RHS[totnode - 1] = vdnpcp[1] # - Place hub pseudo control point freestream influence - # - if length(prescribednodeidxs) == 1 + if iszero(prescribednodeidxs[2]) RHS[totnode + 2] = vdnpcp[2] end diff --git a/src/preprocess/velocities/induced_velocity_matrices.jl b/src/preprocess/velocities/induced_velocity_matrices.jl index 6360cb4c..201686b6 100644 --- a/src/preprocess/velocities/induced_velocity_matrices.jl +++ b/src/preprocess/velocities/induced_velocity_matrices.jl @@ -92,8 +92,8 @@ function induced_velocities_from_vortex_panels_on_points!( enumerate(zip(eachcol(nodemap), influence_length, eachcol(strength))) # Loop through control points being influenced for (i, cpi) in enumerate(eachcol(controlpoint)) - n1 = view(node, :, nmap[1]) - n2 = view(node, :, nmap[2]) + n1 = view(node, :, Int(nmap[1])) + n2 = view(node, :, Int(nmap[2])) # check of self-induced: if isapprox(cpi, 0.5 * (n1 .+ n2)) @@ -126,7 +126,8 @@ function induced_velocities_from_vortex_panels_on_points!( for k in 1:2 # fill the Matrix - VEL[i, nmap[k], :] += gammaj[k] * vel[k, :] + # TODO: is having this gammaj here even necessary, or is it always just 1? If so, really need to just remove it since we're allocating large arrays every time we call this function. + VEL[i, Int(nmap[k]), :] += gammaj[k] * vel[k, :] end #for k end #for i end #for j @@ -222,8 +223,8 @@ function induced_velocities_from_source_panels_on_points!( # Loop through control points being influenced for (i, cpi) in enumerate(eachcol(controlpoint)) - n1 = view(node, :, nmap[1]) - n2 = view(node, :, nmap[2]) + n1 = view(node, :, Int(nmap[1])) + n2 = view(node, :, Int(nmap[2])) # check of self-induced: if isapprox(cpi, 0.5 * (n1 .+ n2)) @@ -256,7 +257,7 @@ function induced_velocities_from_source_panels_on_points!( for k in 1:2 # fill the Matrix - VEL[i, nmap[k], :] += sigmaj[k] * vel[k, :] + VEL[i, Int(nmap[k]), :] += sigmaj[k] * vel[k, :] end #for k end #for i end #for j @@ -360,10 +361,10 @@ function induced_velocities_from_trailing_edge_gap_panel!( for k in 1:2 # fill the Matrix - VEL[i, nmap[k], :] += ndn[k] * vvel[k, :] + ncn[k] * svel[k, :] + VEL[i, Int(nmap[k]), :] += ndn[k] * vvel[k, :] + ncn[k] * svel[k, :] if wake # wake "TE Panels" only have the vortex influence - VEL[i, nmap[k], :] += ndn[k] * vvel[k, :] + VEL[i, Int(nmap[k]), :] += ndn[k] * vvel[k, :] end end #for k end #for j diff --git a/src/process/process.jl b/src/process/process.jl new file mode 100644 index 00000000..55dd00b2 --- /dev/null +++ b/src/process/process.jl @@ -0,0 +1,136 @@ +""" +""" +function process( + solver_options::TS, + solve_parameter_cache_vector, + solve_parameter_cache_dims, + airfoils, + A_bb_LU, + solve_container_caching, + idmaps, + options, +) where {TS<:Union{ExternalSolverOptions,MultiSolverOptions}} + + # - Initialize Aero - # + if options.verbose + println("\nInitializing Velocities") + end + # view the initial conditions out of the inputs cache + solve_parameter_tuple = withdraw_solve_parameter_cache( + solver_options, solve_parameter_cache_vector, solve_parameter_cache_dims + ) + (; vz_rotor, vtheta_rotor, Cm_wake, operating_point, linsys, ivr, ivw) = + solve_parameter_tuple + + # initialize velocities + # TODO; add some sort of unit test for this function + initialize_velocities!( + vz_rotor, + vtheta_rotor, + Cm_wake, + solve_parameter_tuple.operating_point, + (; solve_parameter_tuple.blade_elements..., airfoils...), + (; solve_parameter_tuple.linsys..., A_bb_LU), + solve_parameter_tuple.ivr, + solve_parameter_tuple.ivw, + idmaps.body_totnodes, + idmaps.wake_panel_sheet_be_map, + ) + + # - combine cache and constants - # + const_cache = (; + # - Constants - # + options.verbose, + options.silence_warnings, + options.multipoint_index, + #nlsolve options + solver_options, + # Constant Parameters + airfoils, + A_bb_LU, + idmaps, + solve_parameter_cache_dims, + # Cache(s) + solve_container_caching..., + solve_parameter_tuple..., # contains SIAMFANLEOptions containers needed for solver definition + ) + + # - Solve with ImplicitAD - # + if options.verbose + println("\nSolving Nonlinear System") + end + return ImplicitAD.implicit( + solve, system_residual!, solve_parameter_cache_vector, const_cache + ) +end + +""" +""" +function process( + solver_options::CSORSolverOptions, + solve_parameter_cache_vector, + solve_parameter_cache_dims, + airfoils, + A_bb_LU, + solve_container_caching, + idmaps, + options, +) + + # - Initialize Aero - # + if options.verbose + println("\nInitializing Velocities") + end + # view the initial conditions out of the inputs cache + solve_parameter_tuple = withdraw_solve_parameter_cache( + solver_options, solve_parameter_cache_vector, solve_parameter_cache_dims + ) + (; Gamr, sigr, gamw, operating_point, blade_elements, linsys, ivr, ivw, wakeK) = + solve_parameter_tuple + + # - Initialize States - # + initialize_strengths!( + Gamr, + sigr, + gamw, + operating_point, + (; blade_elements..., airfoils...), + (; linsys..., A_bb_LU), + ivr, + ivw, + wakeK, + idmaps.body_totnodes, + idmaps.wake_nodemap, + idmaps.wake_endnodeidxs, + idmaps.wake_panel_sheet_be_map, + idmaps.wake_node_sheet_be_map, + idmaps.wake_node_ids_along_casing_wake_interface, + idmaps.wake_node_ids_along_centerbody_wake_interface, + ) + + # - combine cache and constants - # + const_cache = (; + # - Constants - # + options.silence_warnings, + options.verbose, + options.multipoint_index, + #CSOR solve options + solver_options, + # Constant Parameters + airfoils, + A_bb_LU, + idmaps, + # Cache(s) + solve_parameter_cache_vector, + solve_parameter_cache_dims, + solve_container_caching..., + ) + + # - Solve with ImplicitAD - # + if options.verbose + println("\nSolving Nonlinear System using CSOR Method") + end + return ImplicitAD.implicit( + solve, CSOR_residual!, solve_parameter_cache_vector, const_cache + ) +end diff --git a/src/process/solve.jl b/src/process/solve.jl index dba7239c..36bee105 100644 --- a/src/process/solve.jl +++ b/src/process/solve.jl @@ -25,6 +25,7 @@ function solve( # General verbose, silence_warnings, + multipoint_index, # nlsolve options solver_options, # Constant Parameters @@ -71,7 +72,7 @@ function solve( TF = eltype(state_variables) resid = MVector{2,TF}(999 * ones(TF, 2)) - conv = solver_options.converged + conv = solver_options.converged[multipoint_index[]] iter = 0 # - SOLVE - # @@ -131,6 +132,7 @@ function solve( # General verbose, silence_warnings, + multipoint_index, # nlsolve options solver_options, # Constant Parameters @@ -208,7 +210,7 @@ function solve( ) # update convergence flag - converged[1] = sol.converged + converged[multipoint_index[]] = sol.converged return sol.minimizer end @@ -225,6 +227,7 @@ function solve( # General verbose, silence_warnings, + multipoint_index, # nlsolve options solver_options, # Constant Parameters @@ -280,7 +283,7 @@ function solve( ) # update convergence flag - solver_options.converged[1] = sol.error <= solver_options.atol + solver_options.converged[multipoint_index[]] = sol.error <= solver_options.atol return sol.x end @@ -301,6 +304,7 @@ function solve( # General verbose, silence_warnings, + multipoint_index, # nlsolve options solver_options, # Constant Parameters @@ -383,7 +387,7 @@ function solve( ) # update convergence flag - converged[1] = result.converged + converged[multipoint_index[]] = result.converged return result.x end @@ -425,7 +429,6 @@ function solve( # - Wrap residual - # if verbose println(" " * "Wrapping Residual") - # println(" " * "Wrapping Residual and Jacobian") end function rwrap!(resid, state_variables, p) @@ -495,7 +498,7 @@ function solve( # update convergence flag # errorcode is 0 if everything is good and could be a bunch of other stuff if not. - solver_options.converged[1] = Bool(iszero(result.errcode)) + solver_options.converged[multipoint_index[]] = Bool(iszero(result.errcode)) return result.solution end @@ -519,6 +522,7 @@ function solve( # General verbose, silence_warnings, + multipoint_index, # nlsolve options solver_options, # Constant Parameters @@ -550,7 +554,6 @@ function solve( # - Wrap residual - # if verbose println(" " * "Wrapping Residual") - # println(" " * "Wrapping Residual and Jacobian") end function rwrap!(resid, state_variables, p) @@ -589,7 +592,7 @@ function solve( ) # update convergence flag - converged[1] = SciMLBase.successful_retcode(sol) + converged[multipoint_index[]] = SciMLBase.successful_retcode(sol) return sol.u end @@ -606,6 +609,7 @@ function solve( # General verbose, silence_warnings, + multipoint_index, # nlsolve options solver_options, # Constant Parameters @@ -681,7 +685,13 @@ function solve( # - SOLVE - # if verbose - println(" " * "Newton Solve Trace:") + if algorithm == :newton + println(" " * "Newton Solve Trace:") + elseif algorithm == :anderson + println(" " * "Anderson Solve Trace:") + else + println(" " * "$(algorithm) Solve Trace:") + end end result = NLsolve.nlsolve( df, @@ -697,7 +707,7 @@ function solve( # update convergence flag # note: need to do this complicated check to ensure that we're only checking |f(x)| 0 - ndp += npanels[end - 1] * 2 + ncp += npanels[end - 1] elseif dte_minus_cbte < 0 ncbp += npanels[end - 1] end + # duct panels are 2x the number of nacelle panels + ndp = 2 * ncp + # duct and center body nodes are 1 more than number of panels ndn = ndp + 1 ncbn = ncbp + 1 @@ -71,6 +74,7 @@ function get_problem_dimensions(paneling_constants) nrotor, # number of rotors nwn, # number of wake nodes nwp, # number of wake panels + ncp, # number of casing panels ndn, # number of duct nodes ncbn, # number of centerbody nodes nbn, # number of body nodes @@ -81,6 +85,7 @@ function get_problem_dimensions(paneling_constants) nwsp, # number of panels in each wake sheet ndwin, # number of duct-wake interfacing nodes ncbwin, # number of centerbody-wake interfacing nodes + nbodies=2, #hard code this for now. ) end @@ -88,14 +93,17 @@ end """ function get_problem_dimensions(body_vortex_panels, rotor_source_panels, wake_vortex_panels) # number of rotors - nrotor = rotor_source_panels.nbodies + nrotor = rotor_source_panels.nbodies[] # number of wake sheets (blade nodes) - nws = wake_vortex_panels.nbodies + nws = wake_vortex_panels.nbodies[] # number of blade elements (panels) nbe = nws - 1 # number of body panels + ncp = + findmin(@view(body_vortex_panels.node[1, 1:Int(body_vortex_panels.nnode[1])]))[2] - + 1 #TODO check this is correct ndp = body_vortex_panels.npanel[1] ncbp = body_vortex_panels.npanel[2] @@ -104,9 +112,9 @@ function get_problem_dimensions(body_vortex_panels, rotor_source_panels, wake_vo ncbn = body_vortex_panels.nnode[2] # number of body panels - nbp = body_vortex_panels.totpanel + nbp = body_vortex_panels.totpanel[] # number of body nodes - nbn = body_vortex_panels.totnode + nbn = body_vortex_panels.totnode[] # number of panels in each wake sheet nwsp = wake_vortex_panels.npanel[1] @@ -114,9 +122,9 @@ function get_problem_dimensions(body_vortex_panels, rotor_source_panels, wake_vo nwsn = wake_vortex_panels.nnode[1] # number of wake panels total - nwp = wake_vortex_panels.totpanel + nwp = wake_vortex_panels.totpanel[] # number of wake nodes total - nwn = wake_vortex_panels.totnode + nwn = wake_vortex_panels.totnode[] # number of duct-wake and centerbody-wake interface nodes ndwin = length( @@ -125,32 +133,32 @@ function get_problem_dimensions(body_vortex_panels, rotor_source_panels, wake_vo in( body_vortex_panels.node[ 1, - body_vortex_panels.endnodeidxs[1, 1]:body_vortex_panels.endnodeidxs[ - 2, 1 - ], + Int(body_vortex_panels.endnodeidxs[1, 1]):Int( + body_vortex_panels.endnodeidxs[2, 1] + ), ], ), wake_vortex_panels.node[ 1, - wake_vortex_panels.endnodeidxs[1, end]:wake_vortex_panels.endnodeidxs[ - 2, end - ], + Int(wake_vortex_panels.endnodeidxs[1, end]):Int( + wake_vortex_panels.endnodeidxs[2, end] + ), ], ), findall( in( body_vortex_panels.node[ 2, - body_vortex_panels.endnodeidxs[1, 1]:body_vortex_panels.endnodeidxs[ - 2, 1 - ], + Int(body_vortex_panels.endnodeidxs[1, 1]):Int( + body_vortex_panels.endnodeidxs[2, 1] + ), ], ), wake_vortex_panels.node[ 2, - wake_vortex_panels.endnodeidxs[1, end]:wake_vortex_panels.endnodeidxs[ - 2, end - ], + Int(wake_vortex_panels.endnodeidxs[1, end]):Int( + wake_vortex_panels.endnodeidxs[2, end] + ), ], ), ), @@ -162,32 +170,32 @@ function get_problem_dimensions(body_vortex_panels, rotor_source_panels, wake_vo in( body_vortex_panels.node[ 1, - body_vortex_panels.endnodeidxs[1, 2]:body_vortex_panels.endnodeidxs[ - 2, 2 - ], + Int(body_vortex_panels.endnodeidxs[1, 2]):Int( + body_vortex_panels.endnodeidxs[2, 2] + ), ], ), wake_vortex_panels.node[ 1, - wake_vortex_panels.endnodeidxs[1, 1]:wake_vortex_panels.endnodeidxs[ - 2, 1 - ], + Int(wake_vortex_panels.endnodeidxs[1, 1]):Int( + wake_vortex_panels.endnodeidxs[2, 1] + ), ], ), findall( in( body_vortex_panels.node[ 2, - body_vortex_panels.endnodeidxs[1, 2]:body_vortex_panels.endnodeidxs[ - 2, 2 - ], + Int(body_vortex_panels.endnodeidxs[1, 2]):Int( + body_vortex_panels.endnodeidxs[2, 2] + ), ], ), wake_vortex_panels.node[ 2, - wake_vortex_panels.endnodeidxs[1, 1]:wake_vortex_panels.endnodeidxs[ - 2, 1 - ], + Int(wake_vortex_panels.endnodeidxs[1, 1]):Int( + wake_vortex_panels.endnodeidxs[2, 1] + ), ], ), ), @@ -197,6 +205,7 @@ function get_problem_dimensions(body_vortex_panels, rotor_source_panels, wake_vo nrotor, # number of rotors nwn, # number of wake nodes nwp, # number of wake panels + ncp, # number of casing panels ndn, # number of duct nodes ncbn, # number of centerbody nodes nbn, # number of body nodes @@ -207,5 +216,6 @@ function get_problem_dimensions(body_vortex_panels, rotor_source_panels, wake_vo nwsp, # number of panels in each wake sheet ndwin, # number of duct-wake interfacing nodes ncbwin, # number of centerbody-wake interfacing nodes + nbodies=2, #hardcode for now. ) end diff --git a/src/utilities/caches.jl b/src/utilities/caches.jl index 496e2f3b..ae049005 100644 --- a/src/utilities/caches.jl +++ b/src/utilities/caches.jl @@ -1,27 +1,263 @@ +""" +""" +function initialize_all_caches(solver_options, paneling_constants) + + # - Pre/Post Containers Cache - # + prepost_container_caching = allocate_prepost_container_cache(paneling_constants) + + # - Solve (Sensitivity) Parameters Cache - # + solve_parameter_caching = allocate_solve_parameter_cache( + solver_options, paneling_constants + ) + + # - Solve Intermediate Containers Cache - # + solve_container_caching = allocate_solve_container_cache( + solver_options, paneling_constants + ) + + return prepost_container_caching, solve_parameter_caching, solve_container_caching +end + ###################################################################### # # # CACHE ALLOCATIONS # # # ###################################################################### +""" +""" +function allocate_body_panel_container(total_length, problem_dimensions) + (; + ndn, # number of duct nodes + ncbn, # number of centerbody nodes + ) = problem_dimensions + + nn = [ndn, ncbn] + # number of panels to generate for each body + np = nn .- 1 + # number of bodies + nb = length(np) + # total number of panels in system + tp = sum(np) + # total number of nodes in system + tn = tp + nb + + return allocate_panel_container(total_length, nn, np, tn, tp, nb) +end + +function allocate_rotor_panel_container(total_length, problem_dimensions) + (; + nrotor, # number of rotors + nws, # number of wake sheets (also rotor nodes) + ) = problem_dimensions + + nn = nws * ones(Int, nrotor) + # number of panels to generate for each body + np = nn .- 1 + # number of bodies + nb = length(np) + # total number of panels in system + tp = sum(np) + # total number of nodes in system + tn = tp + nb + + return allocate_panel_container(total_length, nn, np, tn, tp, nb) +end + +function allocate_wake_panel_container(total_length, problem_dimensions) + (; + nws, # number of wake sheets (also rotor nodes) + nwsn, # number of nodes in each wake sheet + ) = problem_dimensions + + nn = nwsn * ones(Int, nws) + # number of panels to generate for each body + np = nn .- 1 + # number of bodies + nb = length(np) + # total number of panels in system + tp = sum(np) + # total number of nodes in system + tn = tp + nb + + return allocate_panel_container(total_length, nn, np, tn, tp, nb) +end """ """ -function allocate_precomp_container_cache(paneling_constants) - pd = get_problem_dimensions(paneling_constants) +function allocate_panel_container(total_length, nn, np, tn, tp, nb) + s = size(nn) + l = lfs(s) + nnode = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + # number of panels to generate for each body + npanel = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (1,) + l = lfs(s) + nbodies = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (1,) + l = lfs(s) + totpanel = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (1,) + l = lfs(s) + totnode = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (2, tn) + l = lfs(s) + node = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (2, tp) + l = lfs(s) + controlpoint = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + normal = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + tangent = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + nodemap = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (tp) + l = lfs(s) + influence_length = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (nb, 2, 2) + l = lfs(s) + endnodes = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + tenode = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (2, nb) + l = lfs(s) + + itcontrolpoint = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + itnormal = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + ittangent = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + tenormal = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + tendotn = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + tencrossn = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + endnodeidxs = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + endpanelidxs = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + teadjnodeidxs = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (nb,) + l = lfs(s) + teinfluence_length = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + prescribednodeidxs = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + return (; + nnode, + npanel, + nbodies, + totpanel, + totnode, + node, + controlpoint, + normal, + tangent, + nodemap, + influence_length, + endnodes, + tenode, + itcontrolpoint, + itnormal, + ittangent, + tenormal, + tendotn, + tencrossn, + endnodeidxs, + endpanelidxs, + teadjnodeidxs, + teinfluence_length, + prescribednodeidxs, + ), + total_length +end + +""" +""" +function allocate_panel_containers(problem_dimensions, total_length) + body_vortex_panels, total_length = allocate_body_panel_container( + total_length, problem_dimensions + ) + rotor_source_panels, total_length = allocate_rotor_panel_container( + total_length, problem_dimensions + ) + wake_vortex_panels, total_length = allocate_wake_panel_container( + total_length, problem_dimensions + ) + + return (; body_vortex_panels, rotor_source_panels, wake_vortex_panels), total_length +end + +""" +""" +function allocate_prepost_container_cache(paneling_constants::PanelingConstants) + problem_dimensions = get_problem_dimensions(paneling_constants) + + return allocate_prepost_container_cache(problem_dimensions) +end + +function allocate_prepost_container_cache(problem_dimensions) (; + nrotor, # number of rotors + nwn, # number of wake nodes + nwp, # number of wake panels + ncp, # number of casing panels ndn, # number of duct nodes ncbn, # number of centerbody nodes nbn, # number of body nodes nbp, # number of body panels nws, # number of wake sheets (also rotor nodes) - nwsn, # number of nodes per wake sheet - ) = pd + nbe, # number of blade elements (also rotor panels) + nwsn, # number of nodes in each wake sheet + nwsp, # number of panels in each wake sheet + ndwin, # number of duct-wake interfacing nodes + ncbwin, # number of centerbody-wake interfacing nodes + nbodies, + ) = problem_dimensions # - initialize - # total_length = 0 - # - Set up Dimensions - # + ### --- PRE-PROCESSING --- ### + + # - COORDINATES - # dcshape = (2, ndn) dclength = lfs(dcshape) rp_duct_coordinates = (; @@ -41,21 +277,362 @@ function allocate_precomp_container_cache(paneling_constants) wake_grid = (; index=(total_length + 1):(total_length + wglength), shape=wgshape) total_length += wglength - aicnshape = (nbp, nbn) - aicnlength = lfs(aicnshape) - AICn = (; index=(total_length + 1):(total_length + aicnlength), shape=aicnshape) - total_length += aicnlength + rs = (nrotor,) + rl = lfs(rs) + rotor_indices_in_wake = (; index=(total_length + 1):(total_length + rl), shape=rs) + total_length += rl + + # - PANELS - # + panels, total_length = allocate_panel_containers(problem_dimensions, total_length) + + # - INDUCED VELS - # + s = (nbp, nbn, 2) + l = lfs(s) + v_bb = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (nbp, nrotor * nws, 2) + l = lfs(s) + v_br = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (nbp, nwn, 2) + l = lfs(s) + v_bw = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + # - LINSYS - # + s = (nbp, nbn) + l = lfs(s) + AICn = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (2, nbn) + l = lfs(s) + AICpcp = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (nbp,) + l = lfs(s) + vdnb = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (2,) + l = lfs(s) + vdnpcp = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + ##### ----- POST PROCESSING ----- ##### + + ### --- ROTOR Post-Processing Cache --- ### + + s = (nrotor,) + l = lfs(s) + rotor_inviscid_thrust = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_viscous_thrust = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_thrust = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_inviscid_torque = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_viscous_torque = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l - aicpcpshape = (nbp, nbn) - aicpcplength = lfs(aicpcpshape) - AICpcp = (; index=(total_length + 1):(total_length + aicpcplength), shape=aicpcpshape) - total_length += aicpcplength + rotor_torque = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_inviscid_power = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_viscous_power = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_power = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_CT = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_CQ = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_CP = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_efficiency = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + induced_efficiency= (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + + s = (nbe, nrotor) + l = lfs(s) + rotor_inviscid_thrust_dist = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_viscous_thrust_dist = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_inviscid_torque_dist = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_viscous_torque_dist = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_inviscid_power_dist = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + rotor_viscous_power_dist = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + blade_normal_force_per_unit_span = (; + index=(total_length + 1):(total_length + l), shape=s + ) + total_length += l + + blade_tangential_force_per_unit_span = (; + index=(total_length + 1):(total_length + l), shape=s + ) + total_length += l + + cn = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + ct = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + cphi = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + sphi = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + + ### --- BODY Post-Processing Cache --- ### + s = (2, nbp) + l = lfs(s) + Vtot_in = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + Vtot_out = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + Vtot_prejump = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + vtot_body = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + vtot_jump = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + vtot_wake = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + vtot_rotors = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (nbp,) + l = lfs(s) + Vtan_in = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + Vtan_out = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + cp_in = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + cp_out = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (ncp,) + l = lfs(s) + vtan_casing_in = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + vtan_casing_out = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + casing_zpts = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + cp_casing_in = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + cp_casing_out = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (ndn - 1 - ncp,) + l = lfs(s) + vtan_nacelle_in = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + vtan_nacelle_out = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + nacelle_zpts = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + cp_nacelle_in = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + cp_nacelle_out = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (ncbn - 1,) + l = lfs(s) + vtan_centerbody_in = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + vtan_centerbody_out = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + centerbody_zpts = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + cp_centerbody_in = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + cp_centerbody_out = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (ndn - 1,) + l = lfs(s) + duct_jump = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (ncbn - 1,) + l = lfs(s) + centerbody_jump = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (nbp,) + l = lfs(s) + body_jump_term = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + s = (nbodies,) + l = lfs(s) + body_thrust = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + body_force_coefficient = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + ### --- TOTALS Post-Processing Cache --- ### + s = (1,) + l = lfs(s) + total_thrust = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + total_torque = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + total_power = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + total_CT = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + total_CQ = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + total_CP = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + total_efficiency = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l + + ideal_efficiency = (; index=(total_length + 1):(total_length + l), shape=s) + total_length += l # return tuple of initialized cache and associated dimensions return (; - precomp_container_cache=PreallocationTools.DiffCache(zeros(total_length)), - precomp_container_cache_dims=(; - rp_duct_coordinates, rp_centerbody_coordinates, wake_grid, AICn, AICpcp + prepost_container_cache=PreallocationTools.DiffCache(zeros(total_length)), + prepost_container_cache_dims=(; + ### --- PRE --- ### + rp_duct_coordinates, + rp_centerbody_coordinates, + wake_grid, + rotor_indices_in_wake, + panels, + ivb=(; v_bb, v_br, v_bw), + AICn, + AICpcp, + vdnb, + vdnpcp, + ### --- Post --- ### + # - ROTOR - # + rotor_inviscid_thrust, + rotor_inviscid_thrust_dist, + rotor_viscous_thrust, + rotor_viscous_thrust_dist, + rotor_thrust, + rotor_inviscid_torque, + rotor_inviscid_torque_dist, + rotor_viscous_torque, + rotor_viscous_torque_dist, + rotor_torque, + rotor_inviscid_power, + rotor_inviscid_power_dist, + rotor_viscous_power, + rotor_viscous_power_dist, + rotor_power, + rotor_CT, + rotor_CQ, + rotor_CP, + rotor_efficiency, + induced_efficiency, + blade_normal_force_per_unit_span, + blade_tangential_force_per_unit_span, + blade_loading_intermediate_containers=(; cn, ct, cphi, sphi), + # - BODY - # + zpts=(; casing_zpts, nacelle_zpts, centerbody_zpts), + vtan_tuple=(; + Vtot_in, + Vtot_out, + Vtan_in, + Vtan_out, + Vtot_prejump, + vtot_body, + duct_jump, + centerbody_jump, + body_jump_term, + vtot_jump, + vtot_wake, + vtot_rotors, + vtan_casing_in, + vtan_casing_out, + vtan_nacelle_in, + vtan_nacelle_out, + vtan_centerbody_in, + vtan_centerbody_out, + ), + cp_tuple=(; + cp_in, + cp_out, + cp_casing_in, + cp_casing_out, + cp_nacelle_in, + cp_nacelle_out, + cp_centerbody_in, + cp_centerbody_out, + ), + body_thrust, + body_force_coefficient, + # - TOTALS - # + total_thrust, + total_torque, + total_power, + total_CT, + total_CQ, + total_CP, + total_efficiency, + ideal_efficiency, ), ) end @@ -67,13 +644,14 @@ function allocate_solve_parameter_cache( ) # - Get problem dimensions - # - pd = get_problem_dimensions(paneling_constants) + problem_dimensions = get_problem_dimensions(paneling_constants) + return allocate_solve_parameter_cache( - solve_type::CSORSolverOptions, pd; fd_chunk_size=fd_chunk_size, levels=levels + solve_type::CSORSolverOptions, problem_dimensions; fd_chunk_size=fd_chunk_size, levels=levels ) end -function allocate_solve_container_cache_extras(solver_options::SIAMFANLEOptions, input_length, total_length) +function allocate_solve_parameter_extras(solver_options::SIAMFANLEOptions, input_length, total_length) s = (input_length,) l = lfs(s) @@ -93,7 +671,7 @@ function allocate_solve_container_cache_extras(solver_options::SIAMFANLEOptions, return total_length, (; resid_cache_vec, krylov_cache_vec, jvp_cache_vec) end -function allocate_solve_container_cache_extras(solver_options, input_length, total_length) +function allocate_solve_parameter_extras(solver_options, input_length, total_length) return total_length, (;) end @@ -293,16 +871,8 @@ function allocate_solve_parameter_cache( # NOTE: how to know the size of airfoil objects up front? if it's a DFDC type, this is pretty simple, but if it's a CCBlade type, then the size could be anything depending on how many data points are used # airfoils = allocate_airfoil_cache(aftype, nrotor, nbe) - # # - Index Maps - # - # # TODO: decide if you actually need these here. these actually won't/can't change at the optimization level, so you could have a different cache for these - # # TODO: figure out what you renamed these to. Also probably need to just call the function to get the sizes and lengths for some of these - # wake_node_ids_along_casing_wake_interface - # wake_node_ids_along_centerbody_wake_interface - # rotorwakenodeid - # wake_nodemap - # wake_endnodeidxs - # rotor_indices_in_wake - # body_totnodes + # - Index Maps - # + # TODO: would need to figure out how to size these beforehand return (; solve_parameter_cache=PreallocationTools.DiffCache( @@ -343,10 +913,10 @@ function allocate_solve_parameter_cache( ) where {TS<:Union{ExternalSolverOptions, MultiSolverOptions}} # - Get problem dimensions - # - pd = get_problem_dimensions(paneling_constants) + problem_dimensions = get_problem_dimensions(paneling_constants) return allocate_solve_parameter_cache( - solve_type, pd; fd_chunk_size=fd_chunk_size, levels=levels + solve_type, problem_dimensions; fd_chunk_size=fd_chunk_size, levels=levels ) end @@ -362,9 +932,9 @@ function allocate_solve_parameter_cache( ndn, # number of duct nodes ncbn, # number of centerbody nodes nbn, # number of body nodes - nbp, # number of body panels - nws, # number of wake sheets (also rotor nodes) + nbp, # number of body paneallocate_solve_parameter_extrasheets (also rotor nodes) nbe, # number of blade elements (also rotor panels) + nws, #number of wake sheets nwsn, # number of nodes per wake sheet ndwin, # number of duct-wake interfacing nodes ncbwin, # number of centerbody-wake interfacing nodes @@ -387,7 +957,7 @@ function allocate_solve_parameter_cache( Cm_wake = (; index=(total_length + 1):(total_length + l), shape=s) total_length += l -total_length, SIAMFANLE_cache_vecs = allocate_solve_container_cache_extras(solve_type, total_length, total_length) +total_length, SIAMFANLE_cache_vecs = allocate_solve_parameter_extras(solve_type, total_length, total_length) # save state dimensions @@ -570,7 +1140,6 @@ total_length, SIAMFANLE_cache_vecs = allocate_solve_container_cache_extras(solv operating_point=(; Vinf, rhoinf, muinf, asound, Omega), ivr=(; v_rb, v_rr, v_rw), ivw=(; v_wb, v_wr, v_ww), - # ivb=(; v_bb, v_br, v_bw), linsys=(; A_bb, b_bf, A_bw, A_pw, A_br, A_pr), blade_elements=(; B, @@ -596,9 +1165,9 @@ end function allocate_solve_container_cache( solve_type::CSORSolverOptions, paneling_constants::PanelingConstants; fd_chunk_size=12, levels=1 ) - pd = get_problem_dimensions(paneling_constants) + problem_dimensions = get_problem_dimensions(paneling_constants) - return allocate_solve_container_cache(solve_type, pd; fd_chunk_size=fd_chunk_size, levels=levels) + return allocate_solve_container_cache(solve_type, problem_dimensions; fd_chunk_size=fd_chunk_size, levels=levels) end """ @@ -788,10 +1357,10 @@ end function allocate_solve_container_cache( solve_type::TS, paneling_constants::PanelingConstants; fd_chunk_size=12, levels=1 ) where {TS<:Union{ExternalSolverOptions,MultiSolverOptions}} - pd = get_problem_dimensions(paneling_constants) + problem_dimensions = get_problem_dimensions(paneling_constants) return allocate_solve_container_cache( - solve_type, pd; fd_chunk_size=fd_chunk_size, levels=levels + solve_type, problem_dimensions; fd_chunk_size=fd_chunk_size, levels=levels ) end @@ -942,34 +1511,6 @@ function allocate_solve_container_cache( ) end -""" -""" -function allocate_post_cache(paneling_constants) - pd = get_problem_dimensions(paneling_constants) - - (; - nrotor, # number of rotors - nwn, # number of wake nodes - nwp, # number of wake panels - ndn, # number of duct nodes - ncbn, # number of centerbody nodes - nbn, # number of body nodes - nbp, # number of body panels - nws, # number of wake sheets (also rotor nodes) - nbe, # number of blade elements (also rotor panels) - nwsn, # number of nodes per wake sheet - ndwin, # number of duct-wake interfacing nodes - ncbwin, # number of centerbody-wake interfacing nodes - ) = pd - - # - initialize - # - total_length = 0 - - return (; - post_cache=PreallocationTools.DiffCache(zeros(total_length)), post_cache_dims=(;) - ) -end - ###################################################################### # # # CACHE RESHAPING # @@ -978,7 +1519,12 @@ end """ """ -function withdraw_precomp_container_cache(vec, dims) +function withdraw_prepost_container_cache(vec, dims) + + # panels = (;) + # ivb=(; v_bb, v_br, v_bw), + + wake_grid = reshape(@view(vec[dims.wake_grid.index]), dims.wake_grid.shape) rp_duct_coordinates = reshape( @view(vec[dims.rp_duct_coordinates.index]), dims.rp_duct_coordinates.shape ) @@ -986,11 +1532,601 @@ function withdraw_precomp_container_cache(vec, dims) @view(vec[dims.rp_centerbody_coordinates.index]), dims.rp_centerbody_coordinates.shape, ) - wake_grid = reshape(@view(vec[dims.wake_grid.index]), dims.wake_grid.shape) + rotor_indices_in_wake = reshape( + @view(vec[dims.rotor_indices_in_wake.index]), dims.rotor_indices_in_wake.shape + ) AICn = reshape(@view(vec[dims.AICn.index]), dims.AICn.shape) AICpcp = reshape(@view(vec[dims.AICpcp.index]), dims.AICpcp.shape) + vdnb = reshape(@view(vec[dims.vdnb.index]), dims.vdnb.shape) + vdnpcp = reshape(@view(vec[dims.vdnpcp.index]), dims.vdnpcp.shape) - return (; rp_duct_coordinates, rp_centerbody_coordinates, wake_grid, AICn, AICpcp) + ivb = (; + v_bb=reshape(@view(vec[dims.ivb.v_bb.index]), dims.ivb.v_bb.shape), + v_br=reshape(@view(vec[dims.ivb.v_br.index]), dims.ivb.v_br.shape), + v_bw=reshape(@view(vec[dims.ivb.v_bw.index]), dims.ivb.v_bw.shape), + ) + + panels = (; + body_vortex_panels=(; + nnode=reshape( + @view(vec[dims.panels.body_vortex_panels.nnode.index]), + dims.panels.body_vortex_panels.nnode.shape, + ), + npanel=reshape( + @view(vec[dims.panels.body_vortex_panels.npanel.index]), + dims.panels.body_vortex_panels.npanel.shape, + ), + nbodies=reshape( + @view(vec[dims.panels.body_vortex_panels.nbodies.index]), + dims.panels.body_vortex_panels.nbodies.shape, + ), + totpanel=reshape( + @view(vec[dims.panels.body_vortex_panels.totpanel.index]), + dims.panels.body_vortex_panels.totpanel.shape, + ), + totnode=reshape( + @view(vec[dims.panels.body_vortex_panels.totnode.index]), + dims.panels.body_vortex_panels.totnode.shape, + ), + node=reshape( + @view(vec[dims.panels.body_vortex_panels.node.index]), + dims.panels.body_vortex_panels.node.shape, + ), + controlpoint=reshape( + @view(vec[dims.panels.body_vortex_panels.controlpoint.index]), + dims.panels.body_vortex_panels.controlpoint.shape, + ), + normal=reshape( + @view(vec[dims.panels.body_vortex_panels.normal.index]), + dims.panels.body_vortex_panels.normal.shape, + ), + tangent=reshape( + @view(vec[dims.panels.body_vortex_panels.tangent.index]), + dims.panels.body_vortex_panels.tangent.shape, + ), + nodemap=reshape( + @view(vec[dims.panels.body_vortex_panels.nodemap.index]), + dims.panels.body_vortex_panels.nodemap.shape, + ), + influence_length=reshape( + @view(vec[dims.panels.body_vortex_panels.influence_length.index]), + dims.panels.body_vortex_panels.influence_length.shape, + ), + endnodes=reshape( + @view(vec[dims.panels.body_vortex_panels.endnodes.index]), + dims.panels.body_vortex_panels.endnodes.shape, + ), + tenode=reshape( + @view(vec[dims.panels.body_vortex_panels.tenode.index]), + dims.panels.body_vortex_panels.tenode.shape, + ), + itcontrolpoint=reshape( + @view(vec[dims.panels.body_vortex_panels.itcontrolpoint.index]), + dims.panels.body_vortex_panels.itcontrolpoint.shape, + ), + itnormal=reshape( + @view(vec[dims.panels.body_vortex_panels.itnormal.index]), + dims.panels.body_vortex_panels.itnormal.shape, + ), + ittangent=reshape( + @view(vec[dims.panels.body_vortex_panels.ittangent.index]), + dims.panels.body_vortex_panels.ittangent.shape, + ), + tenormal=reshape( + @view(vec[dims.panels.body_vortex_panels.tenormal.index]), + dims.panels.body_vortex_panels.tenormal.shape, + ), + tendotn=reshape( + @view(vec[dims.panels.body_vortex_panels.tendotn.index]), + dims.panels.body_vortex_panels.tendotn.shape, + ), + tencrossn=reshape( + @view(vec[dims.panels.body_vortex_panels.tencrossn.index]), + dims.panels.body_vortex_panels.tencrossn.shape, + ), + endnodeidxs=reshape( + @view(vec[dims.panels.body_vortex_panels.endnodeidxs.index]), + dims.panels.body_vortex_panels.endnodeidxs.shape, + ), + endpanelidxs=reshape( + @view(vec[dims.panels.body_vortex_panels.endpanelidxs.index]), + dims.panels.body_vortex_panels.endpanelidxs.shape, + ), + teadjnodeidxs=reshape( + @view(vec[dims.panels.body_vortex_panels.teadjnodeidxs.index]), + dims.panels.body_vortex_panels.teadjnodeidxs.shape, + ), + teinfluence_length=reshape( + @view(vec[dims.panels.body_vortex_panels.teinfluence_length.index]), + dims.panels.body_vortex_panels.teinfluence_length.shape, + ), + prescribednodeidxs=reshape( + @view(vec[dims.panels.body_vortex_panels.prescribednodeidxs.index]), + dims.panels.body_vortex_panels.prescribednodeidxs.shape, + ), + ), + rotor_source_panels=(; + nnode=reshape( + @view(vec[dims.panels.rotor_source_panels.nnode.index]), + dims.panels.rotor_source_panels.nnode.shape, + ), + npanel=reshape( + @view(vec[dims.panels.rotor_source_panels.npanel.index]), + dims.panels.rotor_source_panels.npanel.shape, + ), + nbodies=reshape( + @view(vec[dims.panels.rotor_source_panels.nbodies.index]), + dims.panels.rotor_source_panels.nbodies.shape, + ), + totpanel=reshape( + @view(vec[dims.panels.rotor_source_panels.totpanel.index]), + dims.panels.rotor_source_panels.totpanel.shape, + ), + totnode=reshape( + @view(vec[dims.panels.rotor_source_panels.totnode.index]), + dims.panels.rotor_source_panels.totnode.shape, + ), + node=reshape( + @view(vec[dims.panels.rotor_source_panels.node.index]), + dims.panels.rotor_source_panels.node.shape, + ), + controlpoint=reshape( + @view(vec[dims.panels.rotor_source_panels.controlpoint.index]), + dims.panels.rotor_source_panels.controlpoint.shape, + ), + normal=reshape( + @view(vec[dims.panels.rotor_source_panels.normal.index]), + dims.panels.rotor_source_panels.normal.shape, + ), + tangent=reshape( + @view(vec[dims.panels.rotor_source_panels.tangent.index]), + dims.panels.rotor_source_panels.tangent.shape, + ), + nodemap=reshape( + @view(vec[dims.panels.rotor_source_panels.nodemap.index]), + dims.panels.rotor_source_panels.nodemap.shape, + ), + influence_length=reshape( + @view(vec[dims.panels.rotor_source_panels.influence_length.index]), + dims.panels.rotor_source_panels.influence_length.shape, + ), + endnodes=reshape( + @view(vec[dims.panels.rotor_source_panels.endnodes.index]), + dims.panels.rotor_source_panels.endnodes.shape, + ), + tenode=reshape( + @view(vec[dims.panels.rotor_source_panels.tenode.index]), + dims.panels.rotor_source_panels.tenode.shape, + ), + itcontrolpoint=reshape( + @view(vec[dims.panels.rotor_source_panels.itcontrolpoint.index]), + dims.panels.rotor_source_panels.itcontrolpoint.shape, + ), + itnormal=reshape( + @view(vec[dims.panels.rotor_source_panels.itnormal.index]), + dims.panels.rotor_source_panels.itnormal.shape, + ), + ittangent=reshape( + @view(vec[dims.panels.rotor_source_panels.ittangent.index]), + dims.panels.rotor_source_panels.ittangent.shape, + ), + tenormal=reshape( + @view(vec[dims.panels.rotor_source_panels.tenormal.index]), + dims.panels.rotor_source_panels.tenormal.shape, + ), + tendotn=reshape( + @view(vec[dims.panels.rotor_source_panels.tendotn.index]), + dims.panels.rotor_source_panels.tendotn.shape, + ), + tencrossn=reshape( + @view(vec[dims.panels.rotor_source_panels.tencrossn.index]), + dims.panels.rotor_source_panels.tencrossn.shape, + ), + endnodeidxs=reshape( + @view(vec[dims.panels.rotor_source_panels.endnodeidxs.index]), + dims.panels.rotor_source_panels.endnodeidxs.shape, + ), + endpanelidxs=reshape( + @view(vec[dims.panels.rotor_source_panels.endpanelidxs.index]), + dims.panels.rotor_source_panels.endpanelidxs.shape, + ), + teadjnodeidxs=reshape( + @view(vec[dims.panels.rotor_source_panels.teadjnodeidxs.index]), + dims.panels.rotor_source_panels.teadjnodeidxs.shape, + ), + teinfluence_length=reshape( + @view(vec[dims.panels.rotor_source_panels.teinfluence_length.index]), + dims.panels.rotor_source_panels.teinfluence_length.shape, + ), + prescribednodeidxs=reshape( + @view(vec[dims.panels.rotor_source_panels.prescribednodeidxs.index]), + dims.panels.rotor_source_panels.prescribednodeidxs.shape, + ), + ), + wake_vortex_panels=(; + nnode=reshape( + @view(vec[dims.panels.wake_vortex_panels.nnode.index]), + dims.panels.wake_vortex_panels.nnode.shape, + ), + npanel=reshape( + @view(vec[dims.panels.wake_vortex_panels.npanel.index]), + dims.panels.wake_vortex_panels.npanel.shape, + ), + nbodies=reshape( + @view(vec[dims.panels.wake_vortex_panels.nbodies.index]), + dims.panels.wake_vortex_panels.nbodies.shape, + ), + totpanel=reshape( + @view(vec[dims.panels.wake_vortex_panels.totpanel.index]), + dims.panels.wake_vortex_panels.totpanel.shape, + ), + totnode=reshape( + @view(vec[dims.panels.wake_vortex_panels.totnode.index]), + dims.panels.wake_vortex_panels.totnode.shape, + ), + node=reshape( + @view(vec[dims.panels.wake_vortex_panels.node.index]), + dims.panels.wake_vortex_panels.node.shape, + ), + controlpoint=reshape( + @view(vec[dims.panels.wake_vortex_panels.controlpoint.index]), + dims.panels.wake_vortex_panels.controlpoint.shape, + ), + normal=reshape( + @view(vec[dims.panels.wake_vortex_panels.normal.index]), + dims.panels.wake_vortex_panels.normal.shape, + ), + tangent=reshape( + @view(vec[dims.panels.wake_vortex_panels.tangent.index]), + dims.panels.wake_vortex_panels.tangent.shape, + ), + nodemap=reshape( + @view(vec[dims.panels.wake_vortex_panels.nodemap.index]), + dims.panels.wake_vortex_panels.nodemap.shape, + ), + influence_length=reshape( + @view(vec[dims.panels.wake_vortex_panels.influence_length.index]), + dims.panels.wake_vortex_panels.influence_length.shape, + ), + endnodes=reshape( + @view(vec[dims.panels.wake_vortex_panels.endnodes.index]), + dims.panels.wake_vortex_panels.endnodes.shape, + ), + tenode=reshape( + @view(vec[dims.panels.wake_vortex_panels.tenode.index]), + dims.panels.wake_vortex_panels.tenode.shape, + ), + itcontrolpoint=reshape( + @view(vec[dims.panels.wake_vortex_panels.itcontrolpoint.index]), + dims.panels.wake_vortex_panels.itcontrolpoint.shape, + ), + itnormal=reshape( + @view(vec[dims.panels.wake_vortex_panels.itnormal.index]), + dims.panels.wake_vortex_panels.itnormal.shape, + ), + ittangent=reshape( + @view(vec[dims.panels.wake_vortex_panels.ittangent.index]), + dims.panels.wake_vortex_panels.ittangent.shape, + ), + tenormal=reshape( + @view(vec[dims.panels.wake_vortex_panels.tenormal.index]), + dims.panels.wake_vortex_panels.tenormal.shape, + ), + tendotn=reshape( + @view(vec[dims.panels.wake_vortex_panels.tendotn.index]), + dims.panels.wake_vortex_panels.tendotn.shape, + ), + tencrossn=reshape( + @view(vec[dims.panels.wake_vortex_panels.tencrossn.index]), + dims.panels.wake_vortex_panels.tencrossn.shape, + ), + endnodeidxs=reshape( + @view(vec[dims.panels.wake_vortex_panels.endnodeidxs.index]), + dims.panels.wake_vortex_panels.endnodeidxs.shape, + ), + endpanelidxs=reshape( + @view(vec[dims.panels.wake_vortex_panels.endpanelidxs.index]), + dims.panels.wake_vortex_panels.endpanelidxs.shape, + ), + teadjnodeidxs=reshape( + @view(vec[dims.panels.wake_vortex_panels.teadjnodeidxs.index]), + dims.panels.wake_vortex_panels.teadjnodeidxs.shape, + ), + teinfluence_length=reshape( + @view(vec[dims.panels.wake_vortex_panels.teinfluence_length.index]), + dims.panels.wake_vortex_panels.teinfluence_length.shape, + ), + prescribednodeidxs=reshape( + @view(vec[dims.panels.wake_vortex_panels.prescribednodeidxs.index]), + dims.panels.wake_vortex_panels.prescribednodeidxs.shape, + ), + ), + ) + + # - ROTOR POST CACHE - # + + rotor_inviscid_thrust = reshape( + @view(vec[dims.rotor_inviscid_thrust.index]), dims.rotor_inviscid_thrust.shape + ) + rotor_inviscid_thrust_dist = reshape( + @view(vec[dims.rotor_inviscid_thrust_dist.index]), + dims.rotor_inviscid_thrust_dist.shape, + ) + rotor_viscous_thrust = reshape( + @view(vec[dims.rotor_viscous_thrust.index]), dims.rotor_viscous_thrust.shape + ) + rotor_viscous_thrust_dist = reshape( + @view(vec[dims.rotor_viscous_thrust_dist.index]), + dims.rotor_viscous_thrust_dist.shape, + ) + rotor_thrust = reshape(@view(vec[dims.rotor_thrust.index]), dims.rotor_thrust.shape) + rotor_inviscid_torque = reshape( + @view(vec[dims.rotor_inviscid_torque.index]), dims.rotor_inviscid_torque.shape + ) + rotor_inviscid_torque_dist = reshape( + @view(vec[dims.rotor_inviscid_torque_dist.index]), + dims.rotor_inviscid_torque_dist.shape, + ) + rotor_viscous_torque = reshape( + @view(vec[dims.rotor_viscous_torque.index]), dims.rotor_viscous_torque.shape + ) + rotor_viscous_torque_dist = reshape( + @view(vec[dims.rotor_viscous_torque_dist.index]), + dims.rotor_viscous_torque_dist.shape, + ) + rotor_torque = reshape(@view(vec[dims.rotor_torque.index]), dims.rotor_torque.shape) + rotor_inviscid_power = reshape( + @view(vec[dims.rotor_inviscid_power.index]), dims.rotor_inviscid_power.shape + ) + rotor_inviscid_power_dist = reshape( + @view(vec[dims.rotor_inviscid_power_dist.index]), + dims.rotor_inviscid_power_dist.shape, + ) + rotor_viscous_power = reshape( + @view(vec[dims.rotor_viscous_power.index]), dims.rotor_viscous_power.shape + ) + rotor_viscous_power_dist = reshape( + @view(vec[dims.rotor_viscous_power_dist.index]), dims.rotor_viscous_power_dist.shape + ) + rotor_power = reshape(@view(vec[dims.rotor_power.index]), dims.rotor_power.shape) + rotor_CT = reshape(@view(vec[dims.rotor_CT.index]), dims.rotor_CT.shape) + rotor_CQ = reshape(@view(vec[dims.rotor_CQ.index]), dims.rotor_CQ.shape) + rotor_CP = reshape(@view(vec[dims.rotor_CP.index]), dims.rotor_CP.shape) + rotor_efficiency = reshape( + @view(vec[dims.rotor_efficiency.index]), dims.rotor_efficiency.shape + ) + induced_efficiency=reshape( + @view(vec[dims.induced_efficiency.index]), + dims.induced_efficiency.shape, + ) + blade_normal_force_per_unit_span = reshape( + @view(vec[dims.blade_normal_force_per_unit_span.index]), + dims.blade_normal_force_per_unit_span.shape, + ) + blade_tangential_force_per_unit_span = reshape( + @view(vec[dims.blade_tangential_force_per_unit_span.index]), + dims.blade_tangential_force_per_unit_span.shape, + ) + blade_loading_intermediate_containers = (; + cn=reshape( + @view(vec[dims.blade_loading_intermediate_containers.cn.index]), + dims.blade_loading_intermediate_containers.cn.shape, + ), + ct=reshape( + @view(vec[dims.blade_loading_intermediate_containers.ct.index]), + dims.blade_loading_intermediate_containers.ct.shape, + ), + cphi=reshape( + @view(vec[dims.blade_loading_intermediate_containers.cphi.index]), + dims.blade_loading_intermediate_containers.cphi.shape, + ), + sphi=reshape( + @view(vec[dims.blade_loading_intermediate_containers.sphi.index]), + dims.blade_loading_intermediate_containers.sphi.shape, + ), + ) + + # - BODY POST CACHE - # + zpts = (; + centerbody_zpts=reshape( + @view(vec[dims.zpts.centerbody_zpts.index]), + dims.zpts.centerbody_zpts.shape, + ), + casing_zpts=reshape( + @view(vec[dims.zpts.casing_zpts.index]), + dims.zpts.casing_zpts.shape, + ), + nacelle_zpts=reshape( + @view(vec[dims.zpts.nacelle_zpts.index]), + dims.zpts.nacelle_zpts.shape, + ), + ) + + vtan_tuple = (; + Vtot_in=reshape( + @view(vec[dims.vtan_tuple.Vtot_in.index]), + dims.vtan_tuple.Vtot_in.shape, + ), + Vtot_out=reshape( + @view(vec[dims.vtan_tuple.Vtot_out.index]), + dims.vtan_tuple.Vtot_out.shape, + ), + Vtan_in=reshape( + @view(vec[dims.vtan_tuple.Vtan_in.index]), + dims.vtan_tuple.Vtan_in.shape, + ), + Vtan_out=reshape( + @view(vec[dims.vtan_tuple.Vtan_out.index]), + dims.vtan_tuple.Vtan_out.shape, + ), + Vtot_prejump=reshape( + @view(vec[dims.vtan_tuple.Vtot_prejump.index]), + dims.vtan_tuple.Vtot_prejump.shape, + ), + vtot_body=reshape( + @view(vec[dims.vtan_tuple.vtot_body.index]), + dims.vtan_tuple.vtot_body.shape, + ), + duct_jump=reshape( + @view(vec[dims.vtan_tuple.duct_jump.index]), + dims.vtan_tuple.duct_jump.shape, + ), + centerbody_jump=reshape( + @view(vec[dims.vtan_tuple.centerbody_jump.index]), + dims.vtan_tuple.centerbody_jump.shape, + ), + body_jump_term=reshape( + @view(vec[dims.vtan_tuple.body_jump_term.index]), + dims.vtan_tuple.body_jump_term.shape, + ), + vtot_jump=reshape( + @view(vec[dims.vtan_tuple.vtot_jump.index]), + dims.vtan_tuple.vtot_jump.shape, + ), + vtot_wake=reshape( + @view(vec[dims.vtan_tuple.vtot_wake.index]), + dims.vtan_tuple.vtot_wake.shape, + ), + vtot_rotors=reshape( + @view(vec[dims.vtan_tuple.vtot_rotors.index]), + dims.vtan_tuple.vtot_rotors.shape, + ), + # Splits: + vtan_casing_in=reshape( + @view(vec[dims.vtan_tuple.vtan_casing_in.index]), + dims.vtan_tuple.vtan_casing_in.shape, + ), + vtan_casing_out=reshape( + @view(vec[dims.vtan_tuple.vtan_casing_out.index]), + dims.vtan_tuple.vtan_casing_out.shape, + ), + vtan_nacelle_in=reshape( + @view(vec[dims.vtan_tuple.vtan_nacelle_in.index]), + dims.vtan_tuple.vtan_nacelle_in.shape, + ), + vtan_nacelle_out=reshape( + @view(vec[dims.vtan_tuple.vtan_nacelle_out.index]), + dims.vtan_tuple.vtan_nacelle_out.shape, + ), + vtan_centerbody_in=reshape( + @view(vec[dims.vtan_tuple.vtan_centerbody_in.index]), + dims.vtan_tuple.vtan_centerbody_in.shape, + ), + vtan_centerbody_out=reshape( + @view(vec[dims.vtan_tuple.vtan_centerbody_out.index]), + dims.vtan_tuple.vtan_centerbody_out.shape, + ), +) + + cp_tuple = (; + cp_in=reshape(@view(vec[dims.cp_tuple.cp_in.index]), dims.cp_tuple.cp_in.shape), + cp_out=reshape(@view(vec[dims.cp_tuple.cp_out.index]), dims.cp_tuple.cp_out.shape), + cp_casing_in=reshape( + @view(vec[dims.cp_tuple.cp_casing_in.index]), dims.cp_tuple.cp_casing_in.shape + ), + cp_casing_out=reshape( + @view(vec[dims.cp_tuple.cp_casing_out.index]), dims.cp_tuple.cp_casing_out.shape + ), + cp_nacelle_in=reshape( + @view(vec[dims.cp_tuple.cp_nacelle_in.index]), dims.cp_tuple.cp_nacelle_in.shape + ), + cp_nacelle_out=reshape( + @view(vec[dims.cp_tuple.cp_nacelle_out.index]), dims.cp_tuple.cp_nacelle_out.shape + ), + cp_centerbody_in=reshape( + @view(vec[dims.cp_tuple.cp_centerbody_in.index]), + dims.cp_tuple.cp_centerbody_in.shape, + ), + cp_centerbody_out=reshape( + @view(vec[dims.cp_tuple.cp_centerbody_out.index]), + dims.cp_tuple.cp_centerbody_out.shape, + ), + ) + body_thrust = reshape( + @view(vec[dims.body_thrust.index]), dims.body_thrust.shape + ) + body_force_coefficient = reshape( + @view(vec[dims.body_force_coefficient.index]), dims.body_force_coefficient.shape + ) + + # - TOTALS POST CACHE - # + total_thrust = reshape( + @view(vec[dims.total_thrust.index]), dims.total_thrust.shape + ) + total_torque = reshape( + @view(vec[dims.total_torque.index]), dims.total_torque.shape + ) + total_power = reshape( + @view(vec[dims.total_power.index]), dims.total_power.shape + ) + total_CT = reshape( + @view(vec[dims.total_CT.index]), dims.total_CT.shape + ) + total_CQ = reshape( + @view(vec[dims.total_CQ.index]), dims.total_CQ.shape + ) + total_CP = reshape( + @view(vec[dims.total_CP.index]), dims.total_CP.shape + ) + total_efficiency = reshape( + @view(vec[dims.total_efficiency.index]), dims.total_efficiency.shape + ) + ideal_efficiency = reshape( + @view(vec[dims.ideal_efficiency.index]), dims.ideal_efficiency.shape + ) + + return (; + # pre + rp_duct_coordinates, + rp_centerbody_coordinates, + wake_grid, + rotor_indices_in_wake, + AICn, + AICpcp, + vdnb, + vdnpcp, + panels, + ivb, + #rotor post + rotor_inviscid_thrust, + rotor_inviscid_thrust_dist, + rotor_viscous_thrust, + rotor_viscous_thrust_dist, + rotor_thrust, + rotor_inviscid_torque, + rotor_inviscid_torque_dist, + rotor_viscous_torque, + rotor_viscous_torque_dist, + rotor_torque, + rotor_inviscid_power, + rotor_inviscid_power_dist, + rotor_viscous_power, + rotor_viscous_power_dist, + rotor_power, + rotor_CT, + rotor_CQ, + rotor_CP, + rotor_efficiency, + induced_efficiency, + blade_normal_force_per_unit_span, + blade_tangential_force_per_unit_span, + blade_loading_intermediate_containers, + # Body Post + zpts, + vtan_tuple, + cp_tuple, + body_thrust, + body_force_coefficient, + # Totals Post + total_thrust, + total_torque, + total_power, + total_CT, + total_CQ, + total_CP, + total_efficiency, + ideal_efficiency, + ) end """ @@ -1035,14 +2171,6 @@ function withdraw_solve_parameter_cache(solver_options::CSORSolverOptions, vec, v_ww=reshape(@view(vec[dims.ivw.v_ww.index]), dims.ivw.v_ww.shape), ) - # # - induced velocities on body - # - # TODO: move to post process withdraw - # ivb = (; - # v_bb=reshape(@view(vec[dims.ivr.v_bb.index]), dims.ivr.v_bb.shape), - # v_br=reshape(@view(vec[dims.ivr.v_br.index]), dims.ivr.v_br.shape), - # v_bw=reshape(@view(vec[dims.ivr.v_bw.index]), dims.ivr.v_bw.shape), - # ) - # - linear system - # linsys = (; A_bb=reshape(@view(vec[dims.linsys.A_bb.index]), dims.linsys.A_bb.shape), @@ -1389,17 +2517,6 @@ function withdraw_solve_container_cache(solver_options::TS, vec, dims) where {TS ) end -""" -""" -function withdraw_post_parameter_cache(vec, dims) - idmaps = (; - ductidsaftofrotors=(), - wake_panel_ids_along_casing_wake_interface=(), - wake_panel_ids_along_centerbody_wake_interface=(), - ) - return (;) -end - #---------------------------------# # Integration Containers # #---------------------------------# diff --git a/src/utilities/types.jl b/src/utilities/types.jl index 5821eeed..6af53fe2 100644 --- a/src/utilities/types.jl +++ b/src/utilities/types.jl @@ -163,7 +163,7 @@ struct GaussLegendre{TN,TW} <: IntegrationMethod weights::TW end -function GaussLegendre(nsamples=20; silence_warnings=true) +function GaussLegendre(nsamples=8; 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 @@ -174,8 +174,8 @@ function GaussLegendre(nsamples=20; silence_warnings=true) end @kwdef struct IntegrationOptions{TN<:IntegrationMethod,TS<:IntegrationMethod} - nominal::TN = GaussLegendre(20) - singular::TS = GaussLegendre(20) + nominal::TN = GaussLegendre(8) + singular::TS = GaussLegendre(8) end #---------------------------------# @@ -298,12 +298,31 @@ end NonlinearSolveOptions(; algorithm=SimpleNonlinearSolve.SimpleNewtonRaphson, atol=1e-12, - additional_kwargs=(; autodiff=SimpleNonlinearSolve.AutoPolyesterForwardDiff()), + additional_kwargs=(; autodiff=SimpleNonlinearSolve.AutoForwardDiff()), ), ] converged::AbstractVector{TB} = [false] end +function ChainSolverOptions(multipoint) + lm = length(multipoint) + return ChainSolverOptions(; + solvers=[ + NLsolveOptions(; algorithm=:anderson, atol=1e-12, converged=fill(false, lm)), + MinpackOptions(; atol=1e-12, converged=fill(false, lm)), + NonlinearSolveOptions(; + algorithm=SimpleNonlinearSolve.SimpleNewtonRaphson, + atol=1e-12, + additional_kwargs=(; + autodiff=SimpleNonlinearSolve.AutoForwardDiff() + ), + converged=fill(false, lm), + ), + ], + converged=fill(false, lm), + ) +end + #---------------------------------# # ELLIPTIC GRID SOLVER TYPES # #---------------------------------# @@ -328,11 +347,21 @@ end # OPTION SET TYPES # #---------------------------------# @kwdef struct Options{ - TB,TF,TS,Tin,TIo<:IntegrationOptions,TSo<:SolverOptionsType,WS<:GridSolverOptionsType + TB, + TBwo, + TF, + TI, + TSf, + TSt, + Tin, + TIo<:IntegrationOptions, + TSo<:SolverOptionsType, + WS<:GridSolverOptionsType, } # - General Options - # verbose::TB = false silence_warnings::TB = true + multipoint_index::TI = [0] # - Geometry Re-interpolation and generation options - # finterp::Tin = FLOWMath.akima autoshiftduct::TB = true @@ -344,20 +373,39 @@ end # - Integration Options - # integration_options::TIo = IntegrationOptions() # - Post-processing Options - # - write_outputs::TB = false - outfile::TS = "outputs.jl" + write_outputs::TBwo = false + outfile::TSf = ["outputs.jl"] checkoutfileexists::TB = false - output_tuple_name::TS = "outs" + output_tuple_name::TSt = "outs" # - Solving Options - # grid_solver_options::WS = GridSolverOptions() - solver_options::TSo = CompositeSolverOptions() + solver_options::TSo = ChainSolverOptions() end function set_options(; kwargs...) return Options(; kwargs...) end -function quicksolve_options(; +function set_options(multipoint; outfile=nothing, output_tuple_name=nothing, kwargs...) + lm = length(multipoint) + + if isnothing(outfile) + outfile = ["outputs" * lpad(i, 3, "0") * ".jl" for i in 1:lm] + end + + if isnothing(output_tuple_name) + output_tuple_name = ["outs" for i in 1:lm] + end + + return set_options(; + solver_options=ChainSolverOptions(multipoint), + outfile=outfile, + output_tuple_name=output_tuple_name, + kwargs..., + ) +end + +function DFDC_options(; grid_solver_options=SLORGridSolverOptions(), solver_options=CSORSolverOptions(), kwargs..., @@ -368,3 +416,28 @@ function quicksolve_options(; kwargs..., ) end + +###################################################################### +# # +# FUNCTION SET HEADER # +# # +###################################################################### + +function promote_propulosor_type(p) + return promote_type( + eltype(p.duct_coordinates), + eltype(p.centerbody_coordinates), + eltype(p.operating_point.Vinf), + eltype(p.operating_point.rhoinf), + eltype(p.operating_point.muinf), + eltype(p.operating_point.asound), + eltype(p.operating_point.Omega), + eltype(p.rotorstator_parameters.B), + eltype(p.rotorstator_parameters.rotorzloc), + eltype(p.rotorstator_parameters.r), + eltype(p.rotorstator_parameters.Rhub), + eltype(p.rotorstator_parameters.Rtip), + eltype(p.rotorstator_parameters.chords), + eltype(p.rotorstator_parameters.twists), + ) +end diff --git a/src/utilities/utils.jl b/src/utilities/utils.jl index 41155709..85f41e84 100644 --- a/src/utilities/utils.jl +++ b/src/utilities/utils.jl @@ -91,26 +91,28 @@ end note: containers must be Arrays, structs of arrays, or tuples of arrays move to utilities """ -function reset_containers!(c) +function reset_containers!(c; exception_keys=[]) if typeof(c) <: AbstractArray #do nothing if it's a string (eltype(c) == String) || (c .= 0) else for p in propertynames(c) - cp = getfield(c, p) - if typeof(cp) <: AbstractArray - if eltype(cp) <: Tuple - for i in 1:length(cp[1]) - for j in 1:length(cp) - cp[j][i] .= 0.0 + if !(p in exception_keys) + cp = getfield(c, p) + if typeof(cp) <: AbstractArray + 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 - #do nothing if it's a string - (eltype(cp) == String) || (cp .= 0) + reset_containers!(cp; exception_keys=exception_keys) end - else - reset_containers!(cp) end end end diff --git a/test/induced_velocities.jl b/test/induced_velocities.jl index 74ef235a..2c9d43cd 100644 --- a/test/induced_velocities.jl +++ b/test/induced_velocities.jl @@ -25,7 +25,7 @@ println("\nINDUCED VELOCITY TESTS") # influencing panel x = [1.0; 2.0] r = [9.0; 10.0] - ip = dt.generate_panels([x'; r']) + ip = dt.generate_panels([x'; r'], isbody=false) #check self influence xi, rho, m, r_influence = dt.calculate_xrm(ip.node, ip.node) @@ -36,7 +36,7 @@ println("\nINDUCED VELOCITY TESTS") #check influence on axis should have influence in x, and zero in r ra = [0.0; 0.0] - ap = dt.generate_panels([x'; ra']) + ap = dt.generate_panels([x'; ra']; isbody=false) xi, rho, m, r_influence = dt.calculate_xrm(ap.controlpoint[:, 1], ip.node[:, 1]) @test dt.vortex_ring_vz(xi, rho, m, r_influence, ip.influence_length[1]) != 0.0 @test dt.vortex_ring_vr(xi, rho, m, r_influence) == 0.0 @@ -52,7 +52,7 @@ println("\nINDUCED VELOCITY TESTS") # influencing panel x = [1.0; 2.0] r = [9.0; 10.0] - ip = dt.generate_panels([x'; r']) + ip = dt.generate_panels([x'; r']; isbody=false) #check self influence xi, rho, m, r_influence = dt.calculate_xrm(ip.node, ip.node) @@ -61,7 +61,7 @@ println("\nINDUCED VELOCITY TESTS") #check influence ON axis should have influence in x, and zero in r ra = [0.0; 0.0] - ap = dt.generate_panels([x'; ra']) + ap = dt.generate_panels([x'; ra']; isbody=false) xi, rho, m, r_influence = dt.calculate_xrm(ap.controlpoint[:, 1], ip.node[:, 1]) @test dt.source_ring_vz(xi, rho, m, r_influence) == 0.0 @test dt.source_ring_vr(xi, rho, m, r_influence) == 0.0 @@ -91,7 +91,7 @@ end gk_integration_options = dt.GaussKronrod() gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) - gl_integration_options = dt.GaussLegendre() + gl_integration_options = dt.GaussLegendre(20) 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) @@ -266,7 +266,7 @@ end gk_integration_options = dt.GaussKronrod() gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) - gl_integration_options = dt.GaussLegendre() + gl_integration_options = dt.GaussLegendre(20) 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) @@ -473,7 +473,7 @@ end gk_integration_options = dt.GaussKronrod() gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) - gl_integration_options = dt.GaussLegendre() + gl_integration_options = dt.GaussLegendre(20) 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) @@ -652,7 +652,7 @@ end gk_integration_options = dt.GaussKronrod() gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) - gl_integration_options = dt.GaussLegendre() + gl_integration_options = dt.GaussLegendre(20) 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) @@ -850,7 +850,7 @@ end gk_integration_options = dt.GaussKronrod() gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) - gl_integration_options = dt.GaussLegendre() + gl_integration_options = dt.GaussLegendre(20) 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) @@ -1045,7 +1045,7 @@ end gk_integration_options = dt.GaussKronrod() gk_cache = dt.allocate_integration_containers(gk_integration_options, 1.0) - gl_integration_options = dt.GaussLegendre() + gl_integration_options = dt.GaussLegendre(20) 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) diff --git a/test/iteration_step_tests.jl b/test/iteration_step_tests.jl index 283708a6..377b159a 100644 --- a/test/iteration_step_tests.jl +++ b/test/iteration_step_tests.jl @@ -21,8 +21,20 @@ println("\nITERATION STEP THROUGH TESTS") include(datapath * "ductape_parameters.jl") #NOTE, for some reason, julia doesn't recognize that this was defined in the above file... operating_point = dt.OperatingPoint(Vinf, rhoinf, muinf, asound, Omega) + propulsor = dt.Propulsor( + rp_duct_coordinates, + rp_centerbody_coordinates, + rotorstator_parameters, + operating_point, + paneling_constants, + reference_parameters, + ) - options = dt.quicksolve_options() + options = dt.DFDC_options(; + integration_options=dt.IntegrationOptions(; + nominal=dt.GaussLegendre(30), singular=dt.GaussLegendre(30) + ), + ) # rotor is at first wake index rotor_indices_in_wake = [1] @@ -42,12 +54,34 @@ println("\nITERATION STEP THROUGH TESTS") )..., ) + # - Set up Pre- and Post-process Cache - # + # Allocate Cache + prepost_container_caching = dt.allocate_prepost_container_cache(problem_dimensions) + + # unpack the caching + (; prepost_container_cache, prepost_container_cache_dims) = prepost_container_caching + + # Get correct cached types + prepost_container_cache_vec = @views PreallocationTools.get_tmp( + prepost_container_cache, Float64(1.0) + ) + + # reset cache + prepost_container_cache_vec .= 0 + + # Reshape Cache + prepost_containers = dt.withdraw_prepost_container_cache( + prepost_container_cache_vec, prepost_container_cache_dims + ) + + # - Set up Solver Sensitivity Paramter Cache - # + # Allocate Cache solve_parameter_caching = dt.allocate_solve_parameter_cache( options.solver_options, problem_dimensions ) - # separate out caching items + # unpack caching (; solve_parameter_cache, solve_parameter_cache_dims) = solve_parameter_caching # get correct cache type @@ -63,12 +97,12 @@ println("\nITERATION STEP THROUGH TESTS") ) # copy over operating point - for f in fieldnames(typeof(operating_point)) - solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) + for f in fieldnames(typeof(propulsor.operating_point)) + solve_parameter_tuple.operating_point[f] .= getfield(propulsor.operating_point, f) end # generate inputs - ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps, panels, problem_dimensions = dt.precompute_parameters!( + ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps, problem_dimensions = dt.precompute_parameters!( solve_parameter_tuple.ivr, solve_parameter_tuple.ivw, solve_parameter_tuple.blade_elements, @@ -82,7 +116,8 @@ println("\nITERATION STEP THROUGH TESTS") paneling_constants, operating_point, reference_parameters, - nothing; #problem dimensions + prepost_containers, + problem_dimensions; grid_solver_options=options.grid_solver_options, autoshiftduct=options.autoshiftduct, itcpshift=options.itcpshift, @@ -111,7 +146,7 @@ println("\nITERATION STEP THROUGH TESTS") ) # - Parameters - # - B = 5 #there are 5 blades + B = [5] #there are 5 blades # get wake sheet length num_wake_z_panels = num_wake_z_nodes - 1 @@ -156,7 +191,7 @@ println("\nITERATION STEP THROUGH TESTS") # get the body strengths from DFDC include(datapath * "iter2_gamb.jl") gamb2 = reformat_sol(res2, nidx) - @test maximum(abs.(solve_containers.gamb .+ gamb2)) < 0.5 + @test maximum(abs.(solve_containers.gamb .+ gamb2)) < 1.5 ##### ----- Test blade element values----- ##### include(datapath * "iter2_blade_element_values.jl") @@ -188,7 +223,7 @@ println("\nITERATION STEP THROUGH TESTS") blade_elements.rotor_panel_centers, ) - @test isapprox(solve_containers.Cz_rotor, be2.Wz, atol=1e-1) + @test isapprox(solve_containers.Cz_rotor, be2.Wz, atol=5e-1) @test isapprox(solve_containers.Cmag_rotor, be2.Wmag, atol=1e-1) @test isapprox(solve_containers.Ctheta_rotor, be2.Wtheta, atol=1e-4) @@ -243,7 +278,7 @@ println("\nITERATION STEP THROUGH TESTS") Gamr2test, deltaG_prev, deltaG, maxBGamr, maxdeltaBGamr, blade_elements.B; test=true ) - @test isapprox(Gamr2test, Gamr2, atol=1e-3) + @test isapprox(Gamr2test, Gamr2, atol=1e-2) ##### ----- Test wake velocities ----- ##### include(datapath * "iter2_wm_wake_panels.jl") @@ -319,7 +354,7 @@ println("\nITERATION STEP THROUGH TESTS") ##### ----- Test Convergence Criteria ----- ##### @test isapprox(maxdeltagamw[], -13.8755674) - @test isapprox(maxdeltaBGamr[], -12.9095955, atol=1e-2) + @test isapprox(maxdeltaBGamr[], -12.9095955, atol=1e-1) @test isapprox(maxBGamr[], 12.3562546) ##### ----- Test ----- ##### diff --git a/test/linear_system_assembly.jl b/test/linear_system_assembly.jl index 42b52a5d..6457a44f 100644 --- a/test/linear_system_assembly.jl +++ b/test/linear_system_assembly.jl @@ -17,7 +17,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") coordinates = [c1, c2] # generate nominal (wake-like) panels - panels = dt.generate_panels(coordinates) + panels = dt.generate_panels(coordinates; isbody=true) # Generate LHS matrix and RHS vector AICn = dt.vortex_aic_boundary_on_boundary( @@ -64,7 +64,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") Vinf = 1.0 #magnitude doesn't matter yet. Vs = Vinf * [1.0 0.0] # axisymmetric, so no radial component - vdnb = dt.freestream_influence_vector(panels.normal, repeat(Vs, panels.totpanel)') + vdnb = dt.freestream_influence_vector(panels.normal, repeat(Vs, panels.totpanel[])') vdnpcp = dt.freestream_influence_vector( panels.itnormal, repeat(Vs, size(panels.itcontrolpoint, 2)) ) @@ -74,8 +74,8 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") AICpcp, panels.npanel, panels.nnode, - panels.totpanel, - panels.totnode, + panels.totpanel[], + panels.totnode[], panels.prescribednodeidxs; dummyval=1.0, ) @@ -84,8 +84,8 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") vdnpcp, panels.npanel, panels.nnode, - panels.totpanel, - panels.totnode, + panels.totpanel[], + panels.totnode[], panels.prescribednodeidxs, ) @@ -148,7 +148,7 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") Vinf = 1.0 #magnitude doesn't matter yet. Vs = Vinf * [1.0 0.0] # axisymmetric, so no radial component - vdnb = dt.freestream_influence_vector(panels.normal, repeat(Vs, panels.totpanel)') + vdnb = dt.freestream_influence_vector(panels.normal, repeat(Vs, panels.totpanel[])') vdnpcp = dt.freestream_influence_vector( panels.itnormal, repeat(Vs, size(panels.itcontrolpoint, 1)) ) @@ -158,8 +158,8 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") AICpcp, panels.npanel, panels.nnode, - panels.totpanel, - panels.totnode, + panels.totpanel[], + panels.totnode[], panels.prescribednodeidxs; dummyval=1.0, ) @@ -168,8 +168,8 @@ println("\nLINEAR SYSTEM ASSEMBLY TESTS") vdnpcp, panels.npanel, panels.nnode, - panels.totpanel, - panels.totnode, + panels.totpanel[], + panels.totnode[], panels.prescribednodeidxs, ) diff --git a/test/panel_generation_tests.jl b/test/panel_generation_tests.jl index c64ebe21..74533e02 100644 --- a/test/panel_generation_tests.jl +++ b/test/panel_generation_tests.jl @@ -10,7 +10,7 @@ println("\nPANEL INITIALIZATION TESTS") # single panel coordinates = [[0.0 1.0; 1.0 1.0]] - panels = dt.generate_panels(coordinates) + panels = dt.generate_panels(coordinates; isbody=false) @test panels.controlpoint[:, 1] == [0.5; 1.0] @test panels.normal[:, 1] == [0.0; 1.0] @@ -23,8 +23,8 @@ println("\nPANEL INITIALIZATION TESTS") @test panels.endnodeidxs == [1; 1;;] @test panels.endpanelidxs == [1; 1;;] @test panels.influence_length == [1.0] - @test panels.totpanel == 1 - @test panels.totnode == 2 + @test panels.totpanel == [1] + @test panels.totnode == [2] @test panels.npanel == [1] @test panels.nnode == [2] @test panels.tenode == [0.0 0.0;;; 0.0 0.0] # this test doesn't make sense in this case, but is present for completeness and to let me know if this behavior of the second index staying zeros changes. @@ -98,8 +98,13 @@ println("\nPANEL INITIALIZATION TESTS") @test panels.tencrossn[:, 1] == [0.0, 0.0] @test all(isapprox.(panels.tencrossn[:, 2], -[sqrt(2) / 2, sqrt(2) / 2])) @test all(isapprox.(panels.teinfluence_length, [0.0, 0.0])) - @test panels.totnode == 8 - @test panels.totpanel == 6 + @test panels.totnode == [8] + @test panels.totpanel == [6] + + ippanels = deepcopy(panels) + dt.reset_containers!(ippanels) + dt.generate_panels!(ippanels, coordinates; isbody=true) + @test all(dt.compare_namedtuples(ippanels,panels)) # - TE panel specific tests - # x1 = [1.0; 0.5; 0.0; 0.5; 1.0] @@ -122,23 +127,23 @@ println("\nPANEL INITIALIZATION TESTS") @test panels.teadjnodeidxs[:, 1] == [5, 1] @test panels.teadjnodeidxs[:, 2] == [8, 8] @test panels.tendotn[1, 1] == - dt.dot(panels.tenormal[:, 1], panels.normal[:, panels.endpanelidxs[2, 1]]) + dt.dot(panels.tenormal[:, 1], panels.normal[:, Int(panels.endpanelidxs[2, 1])]) @test isapprox( panels.tendotn[1, 2], - dt.dot(panels.tenormal[:, 1], panels.normal[:, panels.endpanelidxs[1, 1]]), + dt.dot(panels.tenormal[:, 1], panels.normal[:, Int(panels.endpanelidxs[1, 1])]), ) @test isapprox( panels.tendotn[2, 1], - dt.dot(panels.tenormal[:, 2], panels.normal[:, panels.endpanelidxs[2, 2]]), + dt.dot(panels.tenormal[:, 2], panels.normal[:, Int(panels.endpanelidxs[2, 2])]), ) @test panels.tendotn[2, 2] == 0.0 @test panels.tencrossn[1, 1] == - -dt.cross2mag(panels.tenormal[:, 1], panels.normal[:, panels.endpanelidxs[2, 1]]) + -dt.cross2mag(panels.tenormal[:, 1], panels.normal[:, Int(panels.endpanelidxs[2, 1])]) @test panels.tencrossn[2, 1] == - -dt.cross2mag(panels.tenormal[:, 1], panels.normal[:, panels.endpanelidxs[1, 1]]) + -dt.cross2mag(panels.tenormal[:, 1], panels.normal[:, Int(panels.endpanelidxs[1, 1])]) @test panels.tencrossn[1, 2] == - -dt.cross2mag(panels.tenormal[:, 2], panels.normal[:, panels.endpanelidxs[2, 2]]) + -dt.cross2mag(panels.tenormal[:, 2], panels.normal[:, Int(panels.endpanelidxs[2, 2])]) @test panels.tencrossn[2, 2] == - -dt.cross2mag(panels.tenormal[:, 2], panels.normal[:, panels.endpanelidxs[2, 2]]) + -dt.cross2mag(panels.tenormal[:, 2], panels.normal[:, Int(panels.endpanelidxs[2, 2])]) @test all(isapprox.(panels.teinfluence_length, [0.2, 0.1])) end diff --git a/test/post_processing_tests.jl b/test/post_processing_tests.jl index 407f1596..13fff6c8 100644 --- a/test/post_processing_tests.jl +++ b/test/post_processing_tests.jl @@ -10,20 +10,19 @@ println("\nPOST PROCESSING TESTS") include("data/dfdc_dstestr/DFDC_THRUST_INVISCID.jl") include("data/dfdc_dstestr/DFDC_THRUST_VISCOUS.jl") - B = 5 - Vinf = 0.0 - Vref = 50.0 - rhoinf = 1.226 - Omega = 8000 * pi / 30 + B = [5] + Vinf = [0.0] + Vref = [50.0] + rhoinf = [1.226] + Omega = [8000 * pi / 30] # columns: B*Gamma,Wtheta,dr,section thrust - Wtheta = dfdc_thrust_inviscid[:, 2] + Wtheta = reshape(dfdc_thrust_inviscid[:, 2], (10, 1)) Gamma_tilde = reshape(dfdc_thrust_inviscid[:, 1], (length(Wtheta), 1)) - rotor_panel_length = dfdc_thrust_inviscid[:, 3] - dfdc_tinv_dist = dfdc_thrust_inviscid[:, 4] - - tinv, tinv_dist = dt.inviscid_rotor_trust( - Wtheta, Gamma_tilde, rotor_panel_length, rhoinf + rotor_panel_length = reshape(dfdc_thrust_inviscid[:, 3], (10, 1)) + dfdc_tinv_dist = reshape(dfdc_thrust_inviscid[:, 4], (10, 1)) + tinv, tinv_dist = dt.inviscid_rotor_thrust( + Wtheta, Gamma_tilde, rotor_panel_length, rhoinf[] ) @test all(isapprox.(dfdc_tinv_dist, tinv_dist; atol=1e-6)) @@ -37,7 +36,7 @@ println("\nPOST PROCESSING TESTS") dfdc_tvis_dist = dfdc_thrust_viscous[:, 6] tvis, tvis_dist = dt.viscous_rotor_thrust( - Wx_rotor, Wmag_rotor, B, chord, rotor_panel_length, cd, rhoinf + Wx_rotor, Wmag_rotor, B, chord, rotor_panel_length, cd, rhoinf[] ) @test all(isapprox.(dfdc_tvis_dist, tvis_dist; atol=1e-6)) @@ -51,29 +50,29 @@ println("\nPOST PROCESSING TESTS") include("data/dfdc_dstestr/DFDC_TORQUE_VISCOUS.jl") # columns: B*Gamma,Wx,rpl,rpc,section torque - Wx_rotor = dfdc_torque_inviscid[:, 2] + Wx_rotor = reshape(dfdc_torque_inviscid[:, 2], (10, 1)) Gamma_tilde = reshape(dfdc_torque_inviscid[:, 1], (length(Wx_rotor), 1)) - rpl = dfdc_torque_inviscid[:, 3] - rpc = dfdc_torque_inviscid[:, 4] - dfdc_qinv_dist = dfdc_torque_inviscid[:, 5] + rpl = reshape(dfdc_torque_inviscid[:, 3], (10, 1)) + rpc = reshape(dfdc_torque_inviscid[:, 4], (10, 1)) + dfdc_qinv_dist = reshape(dfdc_torque_inviscid[:, 5], (10, 1)) rotor_inviscid_torque, rotor_inviscid_torque_dist = dt.inviscid_rotor_torque( - Wx_rotor, rpc, rpl, Gamma_tilde, rhoinf + Wx_rotor, rpc, rpl, Gamma_tilde, rhoinf[] ) @test all(isapprox.(dfdc_qinv_dist, rotor_inviscid_torque_dist; atol=1e-6)) # columns: W,chord,cd,Wt,rpl,rpc,section torque - Wmag_rotor = dfdc_torque_viscous[:, 1] - chord = dfdc_torque_viscous[:, 2] - cdrag = dfdc_torque_viscous[:, 3] + Wmag_rotor = reshape(dfdc_torque_viscous[:, 1], (10, 1)) + chord = reshape(dfdc_torque_viscous[:, 2], (10, 1)) + cdrag = reshape(dfdc_torque_viscous[:, 3], (10, 1)) Wtheta_rotor = reshape(dfdc_torque_viscous[:, 4], (length(Wmag_rotor), 1)) - rpl = dfdc_torque_viscous[:, 5] - rpc = dfdc_torque_viscous[:, 6] - dfdc_qvis_dist = dfdc_torque_viscous[:, 7] + rpl = reshape(dfdc_torque_viscous[:, 5], (10, 1)) + rpc = reshape(dfdc_torque_viscous[:, 6], (10, 1)) + dfdc_qvis_dist = reshape(dfdc_torque_viscous[:, 7], (10, 1)) rotor_viscous_torque, rotor_viscous_torque_dist = dt.viscous_rotor_torque( - Wtheta_rotor, Wmag_rotor, B, chord, rpc, rpl, cdrag, rhoinf + Wtheta_rotor, Wmag_rotor, B, chord, rpc, rpl, cdrag, rhoinf[] ) @test all(isapprox.(dfdc_qvis_dist, rotor_viscous_torque_dist; atol=1e-6)) @@ -82,12 +81,12 @@ println("\nPOST PROCESSING TESTS") # ROTOR POWER # #---------------------------------# - pinv = dt.inviscid_rotor_power(rotor_inviscid_torque, Omega) + pinv, _ = dt.rotor_power(rotor_inviscid_torque, rotor_inviscid_torque_dist, Omega) - pvis = dt.viscous_rotor_power(rotor_viscous_torque, Omega) + pvis, _ = dt.rotor_power(rotor_viscous_torque, rotor_viscous_torque_dist, Omega) - @test pvis == rotor_viscous_torque * Omega - @test pinv == rotor_inviscid_torque * Omega + @test pvis[] == (rotor_viscous_torque .* Omega)[] + @test pinv[] == (rotor_inviscid_torque .* Omega)[] #---------------------------------# # DUCT THRUST # @@ -103,14 +102,20 @@ println("\nPOST PROCESSING TESTS") xnormal = dfdc_duct_thrust[:, 4] panel_length = dfdc_duct_thrust[:, 5] panel_radial_position = dfdc_duct_thrust[:, 6] - dfdc_d_thrust = -0.5 * rhoinf * Vref^2 * dfdc_duct_thrust[end, 7] + dfdc_d_thrust = -0.5 * rhoinf[] * Vref[]^2 * dfdc_duct_thrust[end, 7] controlpoint = zeros(2, length(panel_radial_position)) controlpoint[2, :] = panel_radial_position - ductpanel = (; nbodies=1, endpanelidxs=[1; length(panel_length);;], normal=xnormal' .* 1.0, influence_length=panel_length, controlpoint) + ductpanel = (; + nbodies=1, + endpanelidxs=[1; length(panel_length);;], + normal=xnormal' .* 1.0, + influence_length=panel_length, + controlpoint, + ) duct_thrust, _ = dt.forces_from_pressure( - cpleft, cpright, ductpanel; rhoinf=rhoinf, Vref=Vref + cpleft, cpright, ductpanel; rhoinf=rhoinf[], Vref=Vref[] ) @test isapprox(duct_thrust[1], dfdc_d_thrust; atol=1e-4) @@ -121,143 +126,153 @@ println("\nPOST PROCESSING TESTS") xnormal = dfdc_hub_thrust[:, 4] panel_length = dfdc_hub_thrust[:, 5] panel_radial_position = dfdc_hub_thrust[:, 6] - dfdc_h_thrust = -0.5 * rhoinf * Vref^2 * dfdc_hub_thrust[end, 7] + dfdc_h_thrust = -0.5 * rhoinf[] * Vref[]^2 * dfdc_hub_thrust[end, 7] controlpoint = zeros(2, length(panel_radial_position)) controlpoint[2, :] = panel_radial_position - hubpanel = (; nbodies=1, endpanelidxs=[1; length(panel_length);;], normal=xnormal' .* 1.0, influence_length=panel_length, controlpoint) + hubpanel = (; + nbodies=1, + endpanelidxs=[1; length(panel_length);;], + normal=xnormal' .* 1.0, + influence_length=panel_length, + controlpoint, + ) hub_thrust, _ = dt.forces_from_pressure( - cpleft, cpright, hubpanel; rhoinf=rhoinf, Vref=Vref + cpleft, + cpright, + hubpanel; + rhoinf=rhoinf[], + Vref=Vref[], ) @test isapprox(hub_thrust[1], dfdc_h_thrust; atol=1e-4) - -end - -@testset "Body Surface Pressure" begin - datapath = "data/dfdc_init_iter1/" - include(datapath * "ductape_parameters.jl") - include(datapath * "ductape_formatted_dfdc_geometry.jl") - - # generate inputs - inputs = dt.precomputed_inputs( - system_geometry, - rotorstator_parameters, #vector of named tuples - freestream, - reference_parameters; - debug=false, - ) - - # - Parameters - # - B = 5 #there are 5 blades - - # - read in index file - # - # in the case here, element 1 is the hub, 2 is the duct, 3 is the rotor, 4 is the hub wake, 14 is the duct wake, and the rest are the interior wake sheets. - include(datapath * "element_indices.jl") - - # get ranges of element idices - nidx = dfdc_eid[:, 2:3] - pidx = dfdc_eid[:, 4:5] - nidranges = nidx[:, 2] .- nidx[:, 1] .+ 1 - pidranges = pidx[:, 2] .- pidx[:, 1] .+ 1 - - # - read in geometry file to help find indices - # - include(datapath * "dfdc_geometry.jl") - - # get wake sheet length - num_wake_z_panels = num_wake_z_nodes - 1 - - # get number of nodes and panels on bodies interface with wake - nhub_interface_nodes = num_wake_z_nodes - nidranges[4] - nduct_interface_nodes = num_wake_z_nodes - nidranges[end] - nhub_interface_panels = num_wake_z_panels - pidranges[4] - nduct_interface_panels = num_wake_z_panels - pidranges[end] - - # - read in functions to reformat - # - include(datapath * "reformatting_functions.jl") - include(datapath * "iter2_sigr.jl") - include(datapath * "iter1_relaxed_gamw.jl") - gamw1 = reformat_gamw(relaxed_GTH, nidx, nhub_interface_nodes, nduct_interface_nodes) - include(datapath * "iter1_relaxed_bgam.jl") - Gamr1 = reformat_circulation(bgamr1, B) - - # - Body Pressures - # - include(datapath * "iter2_body_cp_withdeltas.jl") - duct_cpR2, hub_cpR2 = reformat_body_cp(body_cpR2, pidx) - - Gamr = zeros(length(Gamr1), 1) .= copy(Gamr1) - sigr = zeros(length(sigr2), 1) .= copy(sigr2) - gamw = -copy(gamw1) - outs = dt.post_process(Gamr, sigr, gamw, inputs; write_outputs=false) - - duct_cp = [outs.bodies.cp_casing_out; outs.bodies.cp_nacelle_out] - centerbody_cp = outs.bodies.cp_centerbody_out - @test all(isapprox.(duct_cp, reverse(duct_cpR2), atol=1e-2)) - @test all(isapprox.(centerbody_cp, reverse(hub_cpR2), atol=1e-2)) end -#TODO: update this test to use DuctAPE.c4b -# @testset "Rotor Loads" begin - -# # set up ccblade example -# Rtip = 10 / 2.0 * 0.0254 # inches to meters -# Rhub = 0.10 * Rtip -# B = 2 # number of blades -# rotor = Rotor(Rhub, Rtip, B) - -# propgeom = [ -# 0.15 0.130 32.76 -# 0.20 0.149 37.19 -# 0.25 0.173 33.54 -# 0.30 0.189 29.25 -# 0.35 0.197 25.64 -# 0.40 0.201 22.54 -# 0.45 0.200 20.27 -# 0.50 0.194 18.46 -# 0.55 0.186 17.05 -# 0.60 0.174 15.97 -# 0.65 0.160 14.87 -# 0.70 0.145 14.09 -# 0.75 0.128 13.39 -# 0.80 0.112 12.84 -# 0.85 0.096 12.25 -# 0.90 0.081 11.37 -# 0.95 0.061 10.19 -# 1.00 0.041 8.99 -# ] - -# r = propgeom[:, 1] * Rtip -# chord = propgeom[:, 2] * Rtip -# theta = propgeom[:, 3] * pi / 180 -# af = AlphaAF("test/data/naca4412.dat") -# sections = Section.(r, chord, theta, Ref(af)) - -# Vinf = 5.0 -# Omega = 5400 * pi / 30 # convert to rad/s -# rho = 1.225 -# op = simple_op.(Vinf, Omega, r, rho) - -# # run ccblade -# out = solve.(Ref(rotor), sections, op) - -# # run ccblade post processing -# T, Q = thrusttorque(rotor, sections, out) -# eff, CT, CQ = nondim(T, Q, Vinf, Omega, rho, rotor, "propeller") - -# # run DuctAPE post processing -# blade_elements = ( -# Omega=Omega, B=B, chords=chord, twists=theta, rbe=r, Rtip=Rtip, Rhub=Rhub -# ) -# aero = dt.get_rotor_loads(out.W, out.phi, out.cl, out.cd, blade_elements, (; Vinf, rho)) -# CTdt = aero.CT -# CQdt = aero.CQ -# effdt = aero.eff -# Npdt = aero.Np -# Tpdt = aero.Tp - -# # compare ccblade and DuctAPE values (should be very similar if not identical) -# @test CT == CTdt -# @test CQ == CQdt -# @test eff == effdt -# end +## TODO: update for new pre-computation functions +#@testset "Body Surface Pressure" begin +# datapath = "data/dfdc_init_iter1/" +# include(datapath * "ductape_parameters.jl") +# include(datapath * "ductape_formatted_dfdc_geometry.jl") + +# # generate inputs +# inputs = dt.precomputed_inputs( +# system_geometry, +# rotorstator_parameters, #vector of named tuples +# freestream, +# reference_parameters; +# debug=false, +# ) + +# # - Parameters - # +# B = [5] #there are 5 blades + +# # - read in index file - # +# # in the case here, element 1 is the hub, 2 is the duct, 3 is the rotor, 4 is the hub wake, 14 is the duct wake, and the rest are the interior wake sheets. +# include(datapath * "element_indices.jl") + +# # get ranges of element idices +# nidx = dfdc_eid[:, 2:3] +# pidx = dfdc_eid[:, 4:5] +# nidranges = nidx[:, 2] .- nidx[:, 1] .+ 1 +# pidranges = pidx[:, 2] .- pidx[:, 1] .+ 1 + +# # - read in geometry file to help find indices - # +# include(datapath * "dfdc_geometry.jl") + +# # get wake sheet length +# num_wake_z_panels = num_wake_z_nodes - 1 + +# # get number of nodes and panels on bodies interface with wake +# nhub_interface_nodes = num_wake_z_nodes - nidranges[4] +# nduct_interface_nodes = num_wake_z_nodes - nidranges[end] +# nhub_interface_panels = num_wake_z_panels - pidranges[4] +# nduct_interface_panels = num_wake_z_panels - pidranges[end] + +# # - read in functions to reformat - # +# include(datapath * "reformatting_functions.jl") +# include(datapath * "iter2_sigr.jl") +# include(datapath * "iter1_relaxed_gamw.jl") +# gamw1 = reformat_gamw(relaxed_GTH, nidx, nhub_interface_nodes, nduct_interface_nodes) +# include(datapath * "iter1_relaxed_bgam.jl") +# Gamr1 = reformat_circulation(bgamr1, B) + +# # - Body Pressures - # +# include(datapath * "iter2_body_cp_withdeltas.jl") +# duct_cpR2, hub_cpR2 = reformat_body_cp(body_cpR2, pidx) + +# Gamr = zeros(length(Gamr1), 1) .= copy(Gamr1) +# sigr = zeros(length(sigr2), 1) .= copy(sigr2) +# gamw = -copy(gamw1) +# outs = dt.post_process(Gamr, sigr, gamw, inputs; write_outputs=false) + +# duct_cp = [outs.bodies.cp_casing_out; outs.bodies.cp_nacelle_out] +# centerbody_cp = outs.bodies.cp_centerbody_out +# @test all(isapprox.(duct_cp, reverse(duct_cpR2), atol=1e-2)) +# @test all(isapprox.(centerbody_cp, reverse(hub_cpR2), atol=1e-2)) +#end + +##TODO: update this test to use DuctAPE.c4b +## @testset "Rotor Loads" begin + +## # set up ccblade example +## Rtip = 10 / 2.0 * 0.0254 # inches to meters +## Rhub = 0.10 * Rtip +## B = 2 # number of blades +## rotor = Rotor(Rhub, Rtip, B) + +## propgeom = [ +## 0.15 0.130 32.76 +## 0.20 0.149 37.19 +## 0.25 0.173 33.54 +## 0.30 0.189 29.25 +## 0.35 0.197 25.64 +## 0.40 0.201 22.54 +## 0.45 0.200 20.27 +## 0.50 0.194 18.46 +## 0.55 0.186 17.05 +## 0.60 0.174 15.97 +## 0.65 0.160 14.87 +## 0.70 0.145 14.09 +## 0.75 0.128 13.39 +## 0.80 0.112 12.84 +## 0.85 0.096 12.25 +## 0.90 0.081 11.37 +## 0.95 0.061 10.19 +## 1.00 0.041 8.99 +## ] + +## r = propgeom[:, 1] * Rtip +## chord = propgeom[:, 2] * Rtip +## theta = propgeom[:, 3] * pi / 180 +## af = AlphaAF("test/data/naca4412.dat") +## sections = Section.(r, chord, theta, Ref(af)) + +## Vinf = 5.0 +## Omega = 5400 * pi / 30 # convert to rad/s +## rho = 1.225 +## op = simple_op.(Vinf, Omega, r, rho) + +## # run ccblade +## out = solve.(Ref(rotor), sections, op) + +## # run ccblade post processing +## T, Q = thrusttorque(rotor, sections, out) +## eff, CT, CQ = nondim(T, Q, Vinf, Omega, rho, rotor, "propeller") + +## # run DuctAPE post processing +## blade_elements = ( +## Omega=Omega, B=B, chords=chord, twists=theta, rbe=r, Rtip=Rtip, Rhub=Rhub +## ) +## aero = dt.get_rotor_loads(out.W, out.phi, out.cl, out.cd, blade_elements, (; Vinf, rho)) +## CTdt = aero.CT +## CQdt = aero.CQ +## effdt = aero.eff +## Npdt = aero.Np +## Tpdt = aero.Tp + +## # compare ccblade and DuctAPE values (should be very similar if not identical) +## @test CT == CTdt +## @test CQ == CQdt +## @test eff == effdt +## end diff --git a/test/pre_processing_tests.jl b/test/pre_processing_tests.jl index 77b8bff2..fdc9d517 100644 --- a/test/pre_processing_tests.jl +++ b/test/pre_processing_tests.jl @@ -353,12 +353,39 @@ end include("data/basic_two_rotor_for_test_NEW.jl") options = dt.set_options() + # - Get Problem Dimensions - # + problem_dimensions = dt.get_problem_dimensions(propulsor.paneling_constants) + + # - Set up Pre- and Post-process Cache - # + # Allocate Cache + prepost_container_caching = dt.allocate_prepost_container_cache( + propulsor.paneling_constants + ) + + # unpack the caching + (; prepost_container_cache, prepost_container_cache_dims) = prepost_container_caching + + # Get correct cached types + prepost_container_cache_vec = @views PreallocationTools.get_tmp( + prepost_container_cache, Float64(1.0) + ) + + # reset cache + prepost_container_cache_vec .= 0 + + # Reshape Cache + prepost_containers = dt.withdraw_prepost_container_cache( + prepost_container_cache_vec, prepost_container_cache_dims + ) + + # - Set up Solver Sensitivity Paramter Cache - # + # Allocate Cache solve_parameter_caching = dt.allocate_solve_parameter_cache( options.solver_options, propulsor.paneling_constants ) - # separate out caching items + # unpack caching (; solve_parameter_cache, solve_parameter_cache_dims) = solve_parameter_caching # get correct cache type @@ -378,15 +405,25 @@ end solve_parameter_tuple.operating_point[f] .= getfield(propulsor.operating_point, f) end - # - Do precomputations - # - ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps, panels, problem_dimensions = dt.precompute_parameters!( + # - Do preprocessutations - # + if options.verbose + println("Pre-computing Parameters") + end + + ##### ----- PERFORM PREPROCESSING COMPUTATIONS ----- ##### + + # - Preprocess - # + ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps, _ = dt.precompute_parameters!( solve_parameter_tuple.ivr, solve_parameter_tuple.ivw, solve_parameter_tuple.blade_elements, solve_parameter_tuple.linsys, solve_parameter_tuple.wakeK, - propulsor; + propulsor, + prepost_containers, + problem_dimensions; grid_solver_options=options.grid_solver_options, + integration_options=options.integration_options, autoshiftduct=options.autoshiftduct, itcpshift=options.itcpshift, axistol=options.axistol, @@ -424,14 +461,41 @@ end # TODO: need to add a better test with more realistic geometry that you can draw more conclusions from # CSOR Solve initialization - options = dt.quicksolve_options() + options = dt.DFDC_options() + + # - Get Problem Dimensions - # + problem_dimensions = dt.get_problem_dimensions(propulsor.paneling_constants) + + # - Set up Pre- and Post-process Cache - # + # Allocate Cache + prepost_container_caching = dt.allocate_prepost_container_cache( + propulsor.paneling_constants + ) + + # unpack the caching + (; prepost_container_cache, prepost_container_cache_dims) = prepost_container_caching + + # Get correct cached types + prepost_container_cache_vec = @views PreallocationTools.get_tmp( + prepost_container_cache, Float64(1.0) + ) + + # reset cache + prepost_container_cache_vec .= 0 + + # Reshape Cache + prepost_containers = dt.withdraw_prepost_container_cache( + prepost_container_cache_vec, prepost_container_cache_dims + ) + + # - Set up Solver Sensitivity Paramter Cache - # # Allocate Cache solve_parameter_caching = dt.allocate_solve_parameter_cache( options.solver_options, propulsor.paneling_constants ) - # separate out caching items + # unpack caching (; solve_parameter_cache, solve_parameter_cache_dims) = solve_parameter_caching # get correct cache type @@ -451,15 +515,25 @@ end solve_parameter_tuple.operating_point[f] .= getfield(propulsor.operating_point, f) end - # - Do precomputations - # - ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps, panels, problem_dimensions = dt.precompute_parameters!( + # - Do preprocessutations - # + if options.verbose + println("Pre-computing Parameters") + end + + ##### ----- PERFORM PREPROCESSING COMPUTATIONS ----- ##### + + # - Preprocess - # + ivb, A_bb_LU, lu_decomp_flag, airfoils, idmaps, _ = dt.precompute_parameters!( solve_parameter_tuple.ivr, solve_parameter_tuple.ivw, solve_parameter_tuple.blade_elements, solve_parameter_tuple.linsys, solve_parameter_tuple.wakeK, - propulsor; + propulsor, + prepost_containers, + problem_dimensions; grid_solver_options=options.grid_solver_options, + integration_options=options.integration_options, autoshiftduct=options.autoshiftduct, itcpshift=options.itcpshift, axistol=options.axistol, diff --git a/test/runtests.jl b/test/runtests.jl index 2f724694..a7a25185 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,6 +40,10 @@ include("state_estimation.jl") include("relaxation_tests.jl") include("wake_aero_tests.jl") +# - post process related tests - # +include("post_processing_tests.jl") + + ########################################################## #EVERYTHING BELOW THIS POINT NEEDS TO BE UPDATED ########################################################## @@ -47,8 +51,6 @@ include("wake_aero_tests.jl") # - Active Development - # # include("solve_checks.jl") -# - post process related tests - # -# include("post_processing_tests.jl") # - Need to update, add, fix, etc. - # # include("body_aero_tests.jl")