Profile Bloch Simulation Method

[rkierulf]

Feature Request

It's possible that the run_spin_precession!() and run_spin_excitation!() functions for the current Bloch simulation method could be made to run faster. For example, inside the main for loop in run_spin_excitation!(), there are some arrays holding intermediate results that could be getting repeatedly allocated:

for s ∈ seq #This iterates over seq, "s = seq[i,:]"
    #Motion
     xt = p.x .+ p.ux(p.x, p.y, p.z, s.t)
     yt = p.y .+ p.uy(p.x, p.y, p.z, s.t)
     zt = p.z .+ p.uz(p.x, p.y, p.z, s.t)
     #Effective field
     ΔBz = p.Δw ./ T(2π * γ) .- s.Δf ./ T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(xt,yt,zt)
     Bz = (s.Gx .* xt .+ s.Gy .* yt .+ s.Gz .* zt) .+ ΔBz
     B = sqrt.(abs.(s.B1) .^ 2 .+ abs.(Bz) .^ 2)
     B[B .== 0] .= eps(T)
     #Spinor Rotation
     φ = T(-2π * γ) * (B .* s.Δt) # TODO: Use trapezoidal integration here (?),  this is just Forward Euler
     mul!( Q(φ, s.B1 ./ B, Bz ./ B), M )
     #Relaxation
     M.xy .= M.xy .* exp.(-s.Δt ./ p.T2)
     M.z  .= M.z  .* exp.(-s.Δt ./ p.T1) .+ p.ρ .* (1 .- exp.(-s.Δt ./ p.T1))
end

It is worth profiling with packages such as BenchmarkingTools.jl (and NVTX.jl for GPU-run simulations) before writing a new kernel-based simulation method.