eth-vaw-glaciology / course-101-0250-00 Goto Github PK

# Lecture 2.2a Parabolic PDEs - solve diffusion equation with explicit forward Euler step
#
# ∂C/∂t = D ∂²C/∂x²

@inbounds @views function diffusion2D_step!(T, Ci, qTx, qTy, dTdt, lam, dt, _dx, _dy)
    qTx     .= .-lam.*(T[2:end  ,2:end-1] .- T[1:end-1,2:end-1]).*_dx   # Fourier's law of heat conduction: qT_x  = -λ ∂T/∂x
    qTy     .= .-lam.*(T[2:end-1,2:end  ] .- T[2:end-1,1:end-1]).*_dy   # ...                               qT_y  = -λ ∂T/∂y
    dTdt    .= Ci[2:end-1,2:end-1].*(                                   # Conservation of energy:           ∂T/∂t = 1/cp (-∂qT_x/∂x - ∂qT_y/∂y)
                 .-(qTx[2:end, :     ] .- qTx[1:end-1, :     ]).*_dx
                 .-(qTy[ :   ,2:end  ] .- qTy[ :     ,1:end-1]).*_dy)
    T[2:end-1,2:end-1] .= T[2:end-1,2:end-1] .+ dt.*dTdt                # Update of temperature             T_new = T_old + ∂t ∂T/∂t
    return nothing
end

@inbounds @views macro d_xa(A) esc(:( ($A[2:end  , :     ] .- $A[1:end-1, :     ]) )) end
@inbounds @views macro d_xi(A) esc(:( ($A[2:end  ,2:end-1] .- $A[1:end-1,2:end-1]) )) end
@inbounds @views macro d_ya(A) esc(:( ($A[ :     ,2:end  ] .- $A[ :     ,1:end-1]) )) end
@inbounds @views macro d_yi(A) esc(:( ($A[2:end-1,2:end  ] .- $A[2:end-1,1:end-1]) )) end
@inbounds @views macro  inn(A) esc(:( $A[2:end-1,2:end-1]                          )) end

@inbounds @views function diffusion2D_step!(T, Ci, qTx, qTy, dTdt, lam, dt, _dx, _dy)
    qTx     .= .-lam.*@d_xi(T).*_dx                              # Fourier's law of heat conduction: qT_x  = -λ ∂T/∂x
    qTy     .= .-lam.*@d_yi(T).*_dy                              # ...                               qT_y  = -λ ∂T/∂y
    dTdt    .= @inn(Ci).*(.-@d_xa(qTx).*_dx .- @d_ya(qTy).*_dy)  # Conservation of energy:           ∂T/∂t = 1/cp (-∂qT_x/∂x - ∂qT_y/∂y)
    @inn(T) .= @inn(T) .+ dt.*dTdt                               # Update of temperature             T_new = T_old + ∂t ∂T/∂t
    return nothing
end