Expose dual tag to user

[EdoAlvarezR]

Description

Expose the tag that is used in dual number with ForwardAD() to the user. This lets the user pass the tag nothing to prevent Julia from compiling the entire Dual chain every time that the objective function is redefined.

If not careful, this could lead to ForwardDiff getting confused with overlapping dual chains (e.g., when the objective function already uses dual numbers for internal computations), so it must be used with caution.

Example

Tested in Julia v1.10.2:

envpath = "./testenvironment"

# -------------------- CREATE THE ENVIRONMENT -----------------------------------
if isdir(envpath)
    rm(envpath, recursive=true)
end
mkpath(envpath)

import Pkg
Pkg.activate(envpath)

Pkg.add("VortexLattice")
Pkg.add(url="https://github.com/byuflowlab/Atmosphere.jl")
Pkg.add("SNOW")


# -------------------- IMPORT MODULES --------------------------------------------
import VortexLattice as vl
import Atmosphere
import SNOW

Pkg.status()


# -------------------- DEFINE AERO MODEL ------------------------------------------

camberlinefun = x -> 0

function evaluate_aero(;
                        magVinf::T1         = 120.0,                         # (m/s) freestream velocity
                        AOA::T2             = 0.0,                           # (deg) angle of attack
                        
                        b::T3               = 10.98061,                      # (m) span length
                        croot::T4           = 0.90445,                       # (m) root chord length
                        taper::T5           = 0.75,                          # taper ratio, ctip/croot
                        
                        xoc_center::T6      = 0.25,                          # x/c center for twist, sweep, and dihedral (0.25 for quarter-chord)
                        dihedral::T7        = 1.0,                           # (deg) dihedral
                        sweep::T8           = 10.0,                          # (deg) sweep
                        
                        twistroot::T9       = 0.0,                           # (deg) twist of root
                        twisttip::T10       = 0.0,                           # (deg) twist of tip

                        camberlinefun       = camberlinefun,                 # Airfoil camberline function
        
                        ns                  = 60,                            # Spanwise elements
                        nc                  = 15,                            # Chordwise elements
                      ) where {T1, T2, T3, T4, T5, T6, T7, T8, T9, T10}

    T = promote_type(T1, T2, T3, T4, T5, T6, T7, T8, T9, T10)


    # Geometry (right half of the wing)
    yle                 = [zero(b), b/2]                        # (m) span stations
    
    chord               = croot .+ croot*(taper-1) * yle/(b/2)  # (m) chord length at each station
    theta               = pi/180 * [twistroot, twisttip]# (rad) twist at each station
    phi                 = pi/180 * dihedral * ones(length(yle)) # (rad) dihedral at each station (I'm not sure what this is for)
    
    xle                 = yle*tand(sweep) - (chord*xoc_center).*cos.(theta) # (m) x-position of TE at each station (x=0 is twist center of root)
    zle                 = yle*tand(dihedral) + (chord*xoc_center).*sin.(theta) # (m) z-position of TE at each station (z=0 is twist center of root)
    
    fc                  = [camberlinefun for i in 1:length(yle)]# Camberline function for each section
    
    
    # discretization parameters
    spacing_s           = vl.Cosine()                   # Spanwise spacing
    spacing_c           = vl.Cosine()                   # Chordwise spacing
    
    
    # reference parameters
    bref = b
    cref = (1 + taper)/2 * croot
    Sref = bref*cref
    rref = [0.0, 0.0, 0.0]                              # (m) reference location for all rotations/moments
    Vinf = magVinf
    ref = vl.Reference(Sref, cref, bref, rref, Vinf)
    
    # freestream parameters
    alpha = pi/180 * AOA                                # (rad) angle of attack
    beta = 0.0                                          # (rad) sideslip angle
    Omega = [0.0; 0.0; 0.0]                             # (rad) vehicle rotation
    fs = vl.Freestream{T}(Vinf, alpha, beta, Omega)
    
    # construct surface
    grid, surface = vl.wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
                                                fc=fc, spacing_s=spacing_s, spacing_c=spacing_c)
    
    # create vector containing all surfaces
    surfaces = [surface]
    
    # perform steady state analysis
    # NOTE: `symmetric=true` assumes the geometry and flow conditions are symmetric about the X-Z axis
    system = vl.steady_analysis(surfaces, ref, fs; symmetric=true)
    
    # retrieve near-field forces
    CF, CM = vl.body_forces(system; frame=vl.Wind())
    
    # perform far-field analysis
    CDiff = vl.far_field_drag(system)
    
    CD, CY, CL = CF
    Cl, Cm, Cn = CM

    return (; CD, CY, CL, Cl, Cm, Cn, Sref)

end



# -------------------- OBJECTIVE AND CONSTRAINTS -----------------------------------------

function obj!(g, X; 
                    # Mission parameters
                    MTOW                = 300 * 0.453592,           # (kg)
                    altitude            = 1000 * 0.3048,             # (m) flight altitude
                    magVinf             = 120.0,                     # (m/s) freestream velocity
                    optargs...
                    )

    # Pre-calculations
    rho, mu, a = Atmosphere.atmospherefit(altitude)

    # Design variables
    variables = (; 
                    twistroot = X[1],
                    twisttip = X[2],
                )

    # Evaluate aero model
    aero = evaluate_aero(; magVinf=magVinf, optargs..., variables...)
    
    # Objective (minimize)
    f = aero.CD / 1e-03

    # Constraints
    CLtrg = MTOW * 9.81 / (0.5*rho*magVinf^2*aero.Sref)     # Target CL

    g[1] = ( aero.CL - CLtrg ) / abs(CLtrg)

    return f
    
end





# -------------------- OPTIMIZATION SETUP -----------------------------------------

# Arguments to aero model
aeroargs = (;
                magVinf             = 120.0,    # (m/s) freestream velocity
                AOA                 = 0.0,      # (deg) angle of attack
    
                b                   = 10.98061, # (m) span length
                croot               = 0.90445,  # (m) root chord length
            )

# Design variables
x0                  = [0.0, -3.0]               # Initial guess
lx                  = [-15.0, -15.0]            # lower bounds on x
ux                  = [15.0, 15.0]              # upper bounds on x

# Constraints
ng                  = 1                         # Number of constraints: 1 == minimum excess CL
lg                  = zeros(ng)                 # Lower bounds on g
ug                  = Inf * ones(ng)            # Upper bounds on g

# Ipopt settings
ip_options = Dict(
                    "max_iter" => 2000,
                    "tol" => 1e-6,
                 )

Running the optimization while we let ForwardAD() automatically define a tag based on the given f function:

# Optimizer options
options             = SNOW.Options( solver=SNOW.IPOPT(ip_options),   # Choosing IPOPT solver
                                    derivatives=SNOW.ForwardAD(),    # AD mode
                                    )

# Wrap objective/constraints function
f = (g, X) -> obj!(g, X; aeroargs...)

# Run optimizer
@time xopt, fopt, info = SNOW.minimize(f, x0, ng, lx, ux, lg, ug, options)

This takes ~140 seconds on my laptop, with most of the time being spent in the first function evaluation since Julia is compiling the objective function with dual numbers tagged to the current in-line definition of f.

Manually setting the dual tag to nothing:

# Optimizer options
options             = SNOW.Options( solver=SNOW.IPOPT(ip_options),   # Choosing IPOPT solver
                                    derivatives=SNOW.ForwardAD(; dualtag=nothing),    # AD mode
                                    )

# Wrap objective/constraints function
f = (g, X) -> obj!(g, X; aeroargs...)

# Run optimizer
@time xopt, fopt, info = SNOW.minimize(f, x0, ng, lx, ux, lg, ug, options)

It still takes ~140 secs the first time we run the optimization, but if we redefine f with new aeroargs... and call SNOW.minimize(f, x0, ng, lx, ux, lg, ug, options) again now it only takes ~8 seconds.