@Elrod: just to confirm that `SIMDDualNumbers.jl`

solved the issue with `ifelse`

in ‘more-than-one’ dimensional function + gradient evaluation with both LV and ForwardDiff.

Thanks again!

```
using SIMDDualNumbers, ForwardDiff, StaticArrays, LoopVectorization, SparseArrays, LinearAlgebra, Plots
function PoissonSparsity!(TC,TW,TE,TS,TN)
TW[2:end-0,:] .= TC[1:end-1,:];
TE[1:end-1,:] .= TC[2:end-0,:];
TS[:,2:end-0] .= TC[:,1:end-1];
TN[:,1:end-1] .= TC[:,2:end-0];
end
function Poisson(x::AbstractVector{<:Real}, data, i, j)
kx, ky, dx, dy, ncx, ncy, TS_BC, TN_BC = data
TC, TW, TE, TS, TN = x[1], x[2], x[3], x[4], x[5]
# Stencil points
TW1 = SIMDDualNumbers.ifelse( i==1, TC, TW )
TE1 = SIMDDualNumbers.ifelse( i==ncx, TC, TE )
TS1 = SIMDDualNumbers.ifelse( j==1, 2TS_BC-TC, TS )
TN1 = SIMDDualNumbers.ifelse( j==ncy, 2TN_BC-TC, TN )
# Fluxes
qxW = -kx*(TC - TW1)/dx
qxE = -kx*(TE1 - TC)/dx
qyS = -ky*(TC - TS1)/dy
qyN = -ky*(TN1 - TC)/dy
# Balance
F = (qxE - qxW)/dx + (qyN - qyS)/dy
return F
end
Poisson(TC::Real,TW::Real,TE::Real,TS::Real,TN::Real,data,i,j) = Poisson(@SVector([TC,TW,TE,TS,TN]),data,i,j) # use multiple dispatch to make it work with the desired argument list
Poisson_grad(Poisson,v,w,x,y,z,d,i,j) = ForwardDiff.gradient(x -> Poisson(x, d, i, j),@SVector([v,w,x,y,z])) # construct gradient function for feta()
function Poisson_valgrad(Poisson,v,w,x,y,z,d,i,j)
sv = @SVector([v,w,x,y,z])
dr = ForwardDiff.gradient!(DiffResults.GradientResult(sv), x -> Poisson(x, d, i, j), sv)
DiffResults.value(dr), DiffResults.gradient(dr)
end
function main()
ncx, ncy = 100, 100
data = (kx=1.0, ky=1.0, dx=0.1, dy=0.1, ncx, ncy, TS_BC=1.0, TN_BC=2.0)
# Arrays
F = fill(0.0, ncx, ncy)
T = fill(0.0, ncx, ncy)
Q = fill(0.0, ncx, ncy)
TC = fill(0.0, ncx, ncy)
TW = fill(0.0, ncx, ncy)
TE = fill(0.0, ncx, ncy)
TS = fill(0.0, ncx, ncy)
TN = fill(0.0, ncx, ncy)
IC = reshape(1:ncx*ncy,ncx,ncy) # equation numbers
IW = fill( 1, ncx, ncy)
IE = fill( 1, ncx, ncy)
IS = fill( 1, ncx, ncy)
IN = fill( 1, ncx, ncy)
# Source term
Q[1:Int(ncx/2),:] .= 1.0
# Compute stencil connectivity
PoissonSparsity!(IC,IW,IE,IS,IN)
# Compute Finite Difference coefficients
PoissonSparsity!(TC.=T,TW,TE,TS,TN)
t1 = @elapsed @turbo for j=1:ncy
for i=1:ncx
F[i,j], (TC[i,j], TW[i,j], TE[i,j], TS[i,j], TN[i,j]) = Poisson_valgrad(Poisson,TC[i,j], TW[i,j], TE[i,j], TS[i,j], TN[i,j],data,i,j)
F[i,j] -= Q[i,j]
end
end
# Residual before solve
println("r = ", norm(F))
# Assemble stiffness matrix
I = fill( 0 , 5*ncx*ncy )
J = fill( 0 , 5*ncx*ncy )
V = fill( 0.0, 5*ncx*ncy )
I .= [IC[:]; IC[:]; IC[:]; IC[:]; IC[:]]
J .= [IC[:]; IW[:]; IE[:]; IS[:]; IN[:]]
V .= [TC[:]; TW[:]; TE[:]; TS[:]; TN[:]]
@time K = sparse(I,J,V)
droptol!(K,1e-6)
# Solve
@time T[:] .-= cholesky(K)\F[:]
# Compute residual after solve
PoissonSparsity!(TC.=T,TW,TE,TS,TN)
t2 = @elapsed @turbo for j=1:ncy
for i=1:ncx
F[i,j] = Poisson(TC[i,j],TW[i,j],TE[i,j],TS[i,j],TN[i,j],data,i,j)
F[i,j] -= Q[i,j]
end
end
println("r = ", norm(F[:]))
# Post-processing
display(plot(heatmap(T',c=:roma)))
println("Time for res. + stencil eval. --> ", t1, " s")
println("Time for res. eval. --> ", t2, " s")
end
for it=1:2
@time main()
end
```