# Implicit steady-state diffusion code using a # residual approach and automatic differentiation using GLMakie using ForwardDiff function compute_q!(q, T, k, Δx) for i in eachindex(q) q[i] = -k[i]*(T[i+1] - T[i])/Δx end return nothing end function residual_T_implicit(Tnew, T, cp, rho, k, Δt, Δx) N = length(T) R = zeros(eltype(Tnew), N) q = zeros(eltype(Tnew), N-1) compute_q!(q, Tnew, k, Δx) for i in 2:length(T)-1 R[i] = cp[i]*rho[i]*(Tnew[i]-T[i])/Δt*0 + (q[i] - q[i-1])/Δx #R[i] = (q[i] - q[i-1])/Δx end # BCs # T[1] = 0.0 # T[end] = 1000.0 R[1 ] = Tnew[1 ] - 0.0 R[end] = Tnew[end] - 1000.0 return R end function main() L = 10000 # length of domain in m n = 100; # number of gridpoints nt = 1000; # number of time steps x = range(0,L, length=n) T = zeros(n) # temperature array Tnew= zeros(n) # new temperature array for PT iterations R = zeros(n) # residual array k = fill(3.0, n-1) # thermal conductivity [W/m/K] #k[1:50] .= 1.5 #k[51:end] .= 3.5 rho = fill(3000, n) # density cp = fill(1000, n) # heat capacity [W/kg/K] Δx = x[2] - x[1] # grid spacing Δt_stable = 0.5*(minimum(rho.*cp) * Δx^2)/ maximum(k) Δt = Δt_stable*1e6 #Δt = 3600*24*365.25*200 T[1] = 0.0 T[end] = 1000.0 Tnew .= T fig,ax,li = lines(x/1e3,T, axis=(xlabel="x [km]", ylabel="T [C]")) # ForwardDiff requires functions of the form R = res(x), or res!(R,x) # It cannot deal with functions like res(x,a,b,c,d), or res!(R, x,a,b,c,d) # So how do we solve this?? # => by using anonymous functions in julia # # res1 = x -> res(x,a,b,c,d) # res2 = (R,x) -> res!(R, x,a,b,c,d) # res1(x) res_T = Tnew -> residual_T_implicit(Tnew, T, cp, rho, k, Δt, Δx) for itime = 1:1 J = ForwardDiff.jacobian(res_T, Tnew) R = res_T(Tnew) ΔT = J\-R Tnew = T + ΔT T .= Tnew if itime % 100 == 0 lines!(ax, x/1e3,T, color=:black) end end lines!(ax, x/1e3,T, color=:black) display(fig) end main()
