Solver Benchmark

The problem of solving the panel strengths that satisfy the no-flow-through condition poses a linear system of equation of the form

\[\begin{align*} A y = b ,\end{align*}\]

where $A$ is the matrix containing the geometry of the panels and wake, $y$ is the vector of panel strengths, and $b$ is the vector of boundary conditions. This is trivially solved as

\[\begin{align*} y = A^{-1}b ,\end{align*}\]

however, depending on the size of $A$ (which depends on the number of panels) it can become inefficient or even unfeasible to explicitely calculate the inverse of $A$. Multiple linear solvers are available in FLOWPanel that avoid explicitely inverting $A$, which are described and benchmarked as follows.

The complete code is available at examples/sweptwing_solverbenchmark.jl but you should also be able to copy and paste these lines after running the first section of this example.

Backslash operator \

Most programming languages implement an operator \ that directly calculates the matrix-vector product $A^{-1}b$. This is more efficient than directly inverting $A$ and then multiplying by $b$, without loosing any accuracy. This linear solver is available under this function:

solve_backslash!(y::AbstractVector, A::AbstractMatrix, b::AbstractVector)

and is used as follows:

import Printf: @sprintf

function calc_lift_drag(body, b, ar, Vinf, magVinf, rho; verbose=true, lbl="")

    CLexp = 0.238
    CDexp = 0.005

    str = ""

    if verbose
        str *= @sprintf "| %15.15s | %-7s | %-7s |\n" "Solver" "CL" "CD"
        str *= "| --------------: | :-----: | :-----: |\n"
    end

    # Calculate velocity away from the body
    Us = pnl.calcfield_U(body, body; fieldname="Uoff",
                            offset=0.02, characteristiclength=(args...)->b/ar)

    # Calculate surface velocity U_∇μ due to the gradient of the doublet strength
    UDeltaGamma = pnl.calcfield_Ugradmu(body)

    # Add both velocities together
    pnl.addfields(body, "Ugradmu", "Uoff")

    # Calculate pressure coeffiecient
    Cps = pnl.calcfield_Cp(body, magVinf; U_fieldname="Uoff")

    # Calculate the force of each panel
    Fs = pnl.calcfield_F(body, magVinf, rho; U_fieldname="Uoff")
    # Calculate total force of the vehicle decomposed as lfit, drag, and sideslip
    Dhat = Vinf/pnl.norm(Vinf)        # Drag direction
    Shat = [0, 1, 0]              # Span direction
    Lhat = pnl.cross(Dhat, Shat)      # Lift direction

    LDS = pnl.calcfield_LDS(body, Lhat, Dhat)

    L = LDS[:, 1]
    D = LDS[:, 2]

    # Force coefficients
    nondim = 0.5*rho*magVinf^2*b^2/ar   # Normalization factor
    CL = sign(pnl.dot(L, Lhat)) * pnl.norm(L) / nondim
    CD = sign(pnl.dot(D, Dhat)) * pnl.norm(D) / nondim
    err = abs(CL-CLexp)/CLexp

    if verbose
        str *= @sprintf "| %15.15s | %-7.4f | %-7.4f |\n" lbl CL CD
        str *= @sprintf "| %15.15s | %-7s | %-7s |\n" "Experimental" CLexp CDexp
        str *= @sprintf "\n\tCL Error:\t%4.3g﹪\n" err*100
    end

    return CL, CD, str

end


# ----- Linear solver test: backslash \ operator
t = @elapsed pnl.solve(body, Uinfs, Das, Dbs; solver=pnl.solve_backslash!)

CL, CD, str = calc_lift_drag(body, b, ar, Vinf, magVinf, rho; lbl="Backslash")

CL Error:	1.89﹪
Run time:	0.44 seconds

LU decomposition

Pre-calculating and re-using the LU decomposition of $A$ is advantageous when the linear system needs to be solved for multiple boundary conditions.

solve_ludiv!(y::AbstractVector,
                A::AbstractMatrix, b::AbstractVector; Alu=nothing)

import FLOWPanel as pnl

Alu = pnl.calc_Alu(A)
pnl.solve_ludiv!(y, A, b; Alu=Alu)

Running the solver:

t = @elapsed pnl.solve(body, Uinfs, Das, Dbs; solver=pnl.solve_ludiv!)

CL, CD, str = calc_lift_drag(body, b, ar, Vinf, magVinf, rho; lbl="LUdiv")

CL Error:	1.89﹪
Run time:	0.44 seconds

Iterative Krylov Solver

Iterative Krylov subspace solvers converge to the right solution rather than directly solving the system of equations. This allows the user to trade off accuracy for computational speed by tuning the tolerance of the solver.

The generalized minimal residual (GMRES) method provided by Krylov.jl is available in FLOWPanel.

solve_gmres!(y::AbstractVector,
                A::AbstractMatrix, b::AbstractVector; Avalue=nothing,
                atol=1e-8, rtol=1e-8, optargs...)

import FLOWPanel as pnl

Avalue = pnl.calc_Avalue(A)
pnl.solve_gmres!(y, A, b; Avalue=Avalue)

Running the solver with tolerance $10^{-8}$:

stats = []                      # Stats of GMRES get stored here

t = @elapsed pnl.solve(body, Uinfs, Das, Dbs;
                        solver = pnl.solve_gmres!,
                        solver_optargs = (atol=1e-8, rtol=1e-8, out=stats))

CL, CD, str = calc_lift_drag(body, b, ar, Vinf, magVinf, rho; lbl="GMRES tol=1e-8")

CL Error:	1.89﹪

Simple stats
 niter: 286
 solved: true
 inconsistent: false
 residuals: []
 Aresiduals: []
 κ₂(A): []
 status: solution good enough given atol and rtol

Run time:	0.71 seconds

Running the solver with tolerance $10^{-2}$:

stats = []                      # Stats of GMRES get stored here

t = @elapsed pnl.solve(body, Uinfs, Das, Dbs;
                        solver = pnl.solve_gmres!,
                        solver_optargs = (atol=1e-2, rtol=1e-2, out=stats))

CL, CD, str = calc_lift_drag(body, b, ar, Vinf, magVinf, rho; lbl="GMRES tol=1e-2")

CL Error:	2.06﹪

Simple stats
 niter: 88
 solved: true
 inconsistent: false
 residuals: []
 Aresiduals: []
 κ₂(A): []
 status: solution good enough given atol and rtol

Run time:	0.39 seconds

Running the solver with tolerance $10^{-1}$:

stats = []                      # Stats of GMRES get stored here

t = @elapsed pnl.solve(body, Uinfs, Das, Dbs;
                        solver = pnl.solve_gmres!,
                        solver_optargs = (atol=1e-1, rtol=1e-1, out=stats))

CL, CD, str = calc_lift_drag(body, b, ar, Vinf, magVinf, rho; lbl="GMRES tol=1e-1")

CL Error:	10.3﹪

Simple stats
 niter: 25
 solved: true
 inconsistent: false
 residuals: []
 Aresiduals: []
 κ₂(A): []
 status: solution good enough given atol and rtol

Run time:	0.29 seconds