Skip GPU excitation relaxation for infinite T1 and T2

KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl

@@ -3,7 +3,7 @@ include("PrecessionKernel.jl")
 include("ExcitationKernel.jl")
 
 """Stores preallocated arrays for use in Bloch GPU run_spin_precession! and run_spin_excitation! functions."""
-struct BlochGPUPrealloc{T} <: PreallocResult{T}
+struct BlochGPUPrealloc{T,SKIP_RELAXATION} <: PreallocResult{T}
     sig_output::AbstractMatrix{Complex{T}}
     sig_output_final::AbstractMatrix{Complex{T}}
     ΔBz::AbstractVector{T}
@@ -18,7 +18,21 @@ function prealloc(
     max_block_length::Integer, 
     groupsize
 ) where {T<:Real, SM<:BlochLikeSimMethods}
-    return BlochGPUPrealloc(
+    return prealloc(
+        sim_method, backend, obj, M, max_block_length, groupsize, Val(false)
+    )
+end
+
+function prealloc(
+    sim_method::SM,
+    backend::KA.GPU,
+    obj::Phantom{T},
+    M::Mag{T},
+    max_block_length::Integer,
+    groupsize,
+    ::Val{SKIP_RELAXATION},
+) where {T<:Real, SM<:BlochLikeSimMethods, SKIP_RELAXATION}
+    return BlochGPUPrealloc{T,SKIP_RELAXATION}(
         KA.zeros(backend, Complex{T}, (cld(size(obj.x, 1), groupsize), max_block_length)),
         KA.zeros(backend, Complex{T}, 1, max_block_length),
         obj.Δw ./ T(2π .* γ)
@@ -33,11 +47,7 @@ prealloc(
     max_block_length::Integer,
     groupsize
 ) where {T<:Real} =
-    BlochGPUPrealloc(
-        KA.zeros(backend, Complex{T}, (cld(size(obj.x, 1), groupsize), max_block_length)),
-        KA.zeros(backend, Complex{T}, 1, max_block_length),
-        obj.Δw ./ T(2π .* γ)
-    )
+    prealloc(sim_method, backend, obj, M, max_block_length, groupsize, Val(false))
 
 prealloc(
     sim_method::BlochMagnusBGL6,
@@ -138,8 +148,8 @@ function run_spin_excitation!(
     sim_method::BlochMagnusBGL4,
     groupsize::Integer,
     backend::KA.Backend,
-    pre::BlochGPUPrealloc
-) where {T<:Real}
+    pre::BlochGPUPrealloc{T,SKIP_RELAXATION}
+) where {T<:Real,SKIP_RELAXATION}
     x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t')
     has_adc = length(sig) > 0
 
@@ -148,7 +158,7 @@ function run_spin_excitation!(
         M.xy, M.z,
         x, y, z, pre.ΔBz, p.T1, p.T2, p.ρ, UInt32(length(M.xy)),
         seq.Gx, seq.Gy, seq.Gz, seq.Δt, seq.Δf, seq.B1, seq.ψ, seq.ADC, UInt32(length(seq.t)),
-        Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc),
+        Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc), Val(SKIP_RELAXATION),
         sim_method,
         ndrange=(cld(length(M.xy), groupsize) * groupsize)
     )
@@ -170,8 +180,8 @@ function run_spin_excitation!(
     sim_method::SM,
     groupsize::Integer,
     backend::KA.Backend,
-    pre::BlochGPUPrealloc
-) where {T<:Real, SM<:BlochLikeSimMethods}
+    pre::BlochGPUPrealloc{T,SKIP_RELAXATION}
+) where {T<:Real, SM<:BlochLikeSimMethods, SKIP_RELAXATION}
     #Motion
     x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t')
     has_adc = length(sig) > 0
@@ -182,7 +192,7 @@ function run_spin_excitation!(
         M.xy, M.z,
         x, y, z, pre.ΔBz, p.T1, p.T2, p.ρ, UInt32(length(M.xy)),
         seq.Gx, seq.Gy, seq.Gz, seq.Δt, seq.Δf, seq.B1, seq.ψ, seq.ADC, UInt32(length(seq.t)),
-        Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc),
+        Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc), Val(SKIP_RELAXATION),
         sim_method,
         ndrange=(cld(length(M.xy), groupsize) * groupsize)
     )
@@ -207,8 +217,8 @@ run_spin_excitation!(
     sim_method::BlochMagnusBGL6,
     groupsize::Integer,
     backend::KA.Backend,
-    pre::BlochGPUPrealloc
-) where {T<:Real} =
+    pre::BlochGPUPrealloc{T,SKIP_RELAXATION}
+) where {T<:Real,SKIP_RELAXATION} =
 begin
     x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t')
     has_adc = length(sig) > 0
@@ -218,7 +228,7 @@ begin
         M.xy, M.z,
         x, y, z, pre.ΔBz, p.T1, p.T2, p.ρ, UInt32(length(M.xy)),
         seq.Gx, seq.Gy, seq.Gz, seq.Δt, seq.Δf, seq.B1, seq.ψ, seq.ADC, UInt32(length(seq.t)),
-        Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc),
+        Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc), Val(SKIP_RELAXATION),
         sim_method,
         ndrange=(cld(length(M.xy), groupsize) * groupsize)
     )

KomaMRICore/src/simulation/SimMethods/Bloch/gpu/ExcitationKernel.jl

@@ -38,9 +38,9 @@ end
     M_xy::AbstractVector{Complex{T}}, M_z, 
     @Const(p_x), @Const(p_y), @Const(p_z), @Const(p_ΔBz), @Const(p_T1), @Const(p_T2), @Const(p_ρ), N_spins,
     @Const(s_Gx), @Const(s_Gy), @Const(s_Gz), @Const(s_Δt), @Const(s_Δf), @Const(s_B1), @Const(s_ψ), @Const(s_ADC), s_length,
-    ::Val{MOTION}, ::Val{USE_WARP_REDUCTION}, ::Val{HAS_ADC},
+    ::Val{MOTION}, ::Val{USE_WARP_REDUCTION}, ::Val{HAS_ADC}, ::Val{SKIP_RELAXATION},
     sim_method::SM
-) where {T, MOTION, USE_WARP_REDUCTION, HAS_ADC, SM <: BlochLikeSimMethods}
+) where {T, MOTION, USE_WARP_REDUCTION, HAS_ADC, SKIP_RELAXATION, SM <: BlochLikeSimMethods}
 
     @uniform N = @groupsize()[1]
     i_l = @index(Local, Linear)
@@ -69,9 +69,11 @@ end
         ΔBz = p_ΔBz[i]
         Mxy_r, Mxy_i = reim(M_xy[i])
         Mz = M_z[i]
-        ρ = p_ρ[i]
-        T1 = p_T1[i]
-        T2 = p_T2[i]
+        if !SKIP_RELAXATION
+            ρ = p_ρ[i]
+            T1 = p_T1[i]
+            T2 = p_T2[i]
+        end
         # Rotating frame -> RF frame
         # M * exp(-i * ψ)
         ψ_start = s_ψ[1]
@@ -103,12 +105,17 @@ end
             Mxy_new_r, Mxy_new_i, Mz_new = rotate_magnetization(θx, θy, θz, Mxy_r, Mxy_i, Mz, T)
             Mxy_new_r, Mxy_new_i, Mz_new = restore_mag_norm(M_norm, Mxy_new_r, Mxy_new_i, Mz_new) # For reduced float precision only.
 
-            # Relaxation
-            E1 = exp(-Δt / T1)
-            E2 = exp(-Δt / T2)
-            Mxy_r = Mxy_new_r * E2
-            Mxy_i = Mxy_new_i * E2
-            Mz = Mz_new * E1 + ρ * (T(1) - E1)
+            if SKIP_RELAXATION
+                Mxy_r = Mxy_new_r
+                Mxy_i = Mxy_new_i
+                Mz = Mz_new
+            else
+                E1 = exp(-Δt / T1)
+                E2 = exp(-Δt / T2)
+                Mxy_r = Mxy_new_r * E2
+                Mxy_i = Mxy_new_i * E2
+                Mz = Mz_new * E1 + ρ * (T(1) - E1)
+            end
 
             Bx_0, By_0, Bz_0 = Bx_1, By_1, Bz_1
         end

KomaMRICore/src/simulation/SimMethods/BlochMagnus/gpu/MagnusBGLKernel.jl

@@ -5,9 +5,9 @@
     M_xy::AbstractVector{Complex{T}}, M_z,
     @Const(p_x), @Const(p_y), @Const(p_z), @Const(p_ΔBz), @Const(p_T1), @Const(p_T2), @Const(p_ρ), N_spins,
     @Const(s_Gx), @Const(s_Gy), @Const(s_Gz), @Const(s_Δt), @Const(s_Δf), @Const(s_B1), @Const(s_ψ), @Const(s_ADC), s_length,
-    ::Val{MOTION}, ::Val{USE_WARP_REDUCTION}, ::Val{HAS_ADC},
+    ::Val{MOTION}, ::Val{USE_WARP_REDUCTION}, ::Val{HAS_ADC}, ::Val{SKIP_RELAXATION},
     sim_method::SM
-) where {T, MOTION, USE_WARP_REDUCTION, HAS_ADC, SM<:Union{BlochMagnusBGL4,BlochMagnusBGL6}}
+) where {T, MOTION, USE_WARP_REDUCTION, HAS_ADC, SKIP_RELAXATION, SM<:Union{BlochMagnusBGL4,BlochMagnusBGL6}}
 
     @uniform N = @groupsize()[1]
     i_l = @index(Local, Linear)
@@ -33,9 +33,11 @@
         ΔBz = p_ΔBz[i]
         Mxy_r, Mxy_i = reim(M_xy[i])
         Mz = M_z[i]
-        ρ = p_ρ[i]
-        T1 = p_T1[i]
-        T2 = p_T2[i]
+        if !SKIP_RELAXATION
+            ρ = p_ρ[i]
+            T1 = p_T1[i]
+            T2 = p_T2[i]
+        end
 
         ψ_start = s_ψ[1]
         if !iszero(ψ_start)
@@ -75,14 +77,16 @@
                 sim_method,
             )
             M_norm = mag_norm(T, Mxy_r, Mxy_i, Mz)
-            Mxy_new_r, Mxy_new_i, Mz_new = rotate_magnetization(θx, θy, θz, Mxy_r, Mxy_i, Mz, T)
-            Mxy_new_r, Mxy_new_i, Mz_new = restore_mag_norm(M_norm, Mxy_new_r, Mxy_new_i, Mz_new) # For reduced float precision only.
-
-            E1 = exp(-Δt / T1)
-            E2 = exp(-Δt / T2)
-            Mxy_r = Mxy_new_r * E2
-            Mxy_i = Mxy_new_i * E2
-            Mz = Mz_new * E1 + ρ * (T(1) - E1)
+            Mxy_r, Mxy_i, Mz = rotate_magnetization(θx, θy, θz, Mxy_r, Mxy_i, Mz, T)
+            Mxy_r, Mxy_i, Mz = restore_mag_norm(M_norm, Mxy_r, Mxy_i, Mz) # For reduced float precision only.
+
+            if !SKIP_RELAXATION
+                E1 = exp(-Δt / T1)
+                E2 = exp(-Δt / T2)
+                Mxy_r *= E2
+                Mxy_i *= E2
+                Mz = Mz * E1 + ρ * (T(1) - E1)
+            end
         end
 
         if HAS_ADC && s_ADC[s_end]

KomaMRICore/src/simulation/SimMethods/BlochMagnus/gpu/MagnusGLKernel.jl

@@ -5,9 +5,9 @@
     M_xy::AbstractVector{Complex{T}}, M_z,
     @Const(p_x), @Const(p_y), @Const(p_z), @Const(p_ΔBz), @Const(p_T1), @Const(p_T2), @Const(p_ρ), N_spins,
     @Const(s_Gx), @Const(s_Gy), @Const(s_Gz), @Const(s_Δt), @Const(s_Δf), @Const(s_B1), @Const(s_ψ), @Const(s_ADC), s_length,
-    ::Val{MOTION}, ::Val{USE_WARP_REDUCTION}, ::Val{HAS_ADC},
+    ::Val{MOTION}, ::Val{USE_WARP_REDUCTION}, ::Val{HAS_ADC}, ::Val{SKIP_RELAXATION},
     sim_method::SM
-) where {T, MOTION, USE_WARP_REDUCTION, HAS_ADC, SM<:Union{BlochMagnusGL2,BlochMagnusGL4}}
+) where {T, MOTION, USE_WARP_REDUCTION, HAS_ADC, SKIP_RELAXATION, SM<:Union{BlochMagnusGL2,BlochMagnusGL4}}
 
     @uniform N = @groupsize()[1]
     i_l = @index(Local, Linear)
@@ -33,9 +33,11 @@
         ΔBz = p_ΔBz[i]
         Mxy_r, Mxy_i = reim(M_xy[i])
         Mz = M_z[i]
-        ρ = p_ρ[i]
-        T1 = p_T1[i]
-        T2 = p_T2[i]
+        if !SKIP_RELAXATION
+            ρ = p_ρ[i]
+            T1 = p_T1[i]
+            T2 = p_T2[i]
+        end
 
         ψ_start = s_ψ[1]
         if !iszero(ψ_start)
@@ -70,14 +72,16 @@
                 sim_method,
             )
             M_norm = mag_norm(T, Mxy_r, Mxy_i, Mz)
-            Mxy_new_r, Mxy_new_i, Mz_new = rotate_magnetization(θx, θy, θz, Mxy_r, Mxy_i, Mz, T)
-            Mxy_new_r, Mxy_new_i, Mz_new = restore_mag_norm(M_norm, Mxy_new_r, Mxy_new_i, Mz_new) # For reduced float precision only.
-
-            E1 = exp(-Δt / T1)
-            E2 = exp(-Δt / T2)
-            Mxy_r = Mxy_new_r * E2
-            Mxy_i = Mxy_new_i * E2
-            Mz = Mz_new * E1 + ρ * (T(1) - E1)
+            Mxy_r, Mxy_i, Mz = rotate_magnetization(θx, θy, θz, Mxy_r, Mxy_i, Mz, T)
+            Mxy_r, Mxy_i, Mz = restore_mag_norm(M_norm, Mxy_r, Mxy_i, Mz) # For reduced float precision only.
+
+            if !SKIP_RELAXATION
+                E1 = exp(-Δt / T1)
+                E2 = exp(-Δt / T2)
+                Mxy_r *= E2
+                Mxy_i *= E2
+                Mz = Mz * E1 + ρ * (T(1) - E1)
+            end
         end
 
         if HAS_ADC && s_ADC[s_end]

KomaMRICore/src/simulation/SimMethods/BlochMagnus/gpu/MagnusMidKernel.jl

@@ -5,9 +5,9 @@
     M_xy::AbstractVector{Complex{T}}, M_z,
     @Const(p_x), @Const(p_y), @Const(p_z), @Const(p_ΔBz), @Const(p_T1), @Const(p_T2), @Const(p_ρ), N_spins,
     @Const(s_Gx), @Const(s_Gy), @Const(s_Gz), @Const(s_Δt), @Const(s_Δf), @Const(s_B1), @Const(s_ψ), @Const(s_ADC), s_length,
-    ::Val{MOTION}, ::Val{USE_WARP_REDUCTION}, ::Val{HAS_ADC},
+    ::Val{MOTION}, ::Val{USE_WARP_REDUCTION}, ::Val{HAS_ADC}, ::Val{SKIP_RELAXATION},
     sim_method::BlochMagnusMid2
-) where {T, MOTION, USE_WARP_REDUCTION, HAS_ADC}
+) where {T, MOTION, USE_WARP_REDUCTION, HAS_ADC, SKIP_RELAXATION}
 
     @uniform N = @groupsize()[1]
     i_l = @index(Local, Linear)
@@ -33,9 +33,11 @@
         ΔBz = p_ΔBz[i]
         Mxy_r, Mxy_i = reim(M_xy[i])
         Mz = M_z[i]
-        ρ = p_ρ[i]
-        T1 = p_T1[i]
-        T2 = p_T2[i]
+        if !SKIP_RELAXATION
+            ρ = p_ρ[i]
+            T1 = p_T1[i]
+            T2 = p_T2[i]
+        end
 
         ψ_start = s_ψ[1]
         if !iszero(ψ_start)
@@ -60,14 +62,16 @@
             Δt = s_Δt[s_idx] + s_Δt[s_mid]
             θx, θy, θz = rotation_vector(Bx_m, By_m, Bz_m, Δt, sim_method)
             M_norm = mag_norm(T, Mxy_r, Mxy_i, Mz)
-            Mxy_new_r, Mxy_new_i, Mz_new = rotate_magnetization(θx, θy, θz, Mxy_r, Mxy_i, Mz, T)
-            Mxy_new_r, Mxy_new_i, Mz_new = restore_mag_norm(M_norm, Mxy_new_r, Mxy_new_i, Mz_new) # For reduced float precision only.
+            Mxy_r, Mxy_i, Mz = rotate_magnetization(θx, θy, θz, Mxy_r, Mxy_i, Mz, T)
+            Mxy_r, Mxy_i, Mz = restore_mag_norm(M_norm, Mxy_r, Mxy_i, Mz) # For reduced float precision only.
 
-            E1 = exp(-Δt / T1)
-            E2 = exp(-Δt / T2)
-            Mxy_r = Mxy_new_r * E2
-            Mxy_i = Mxy_new_i * E2
-            Mz = Mz_new * E1 + ρ * (T(1) - E1)
+            if !SKIP_RELAXATION
+                E1 = exp(-Δt / T1)
+                E2 = exp(-Δt / T2)
+                Mxy_r *= E2
+                Mxy_i *= E2
+                Mz = Mz * E1 + ρ * (T(1) - E1)
+            end
         end
 
         if HAS_ADC && s_ADC[s_end]

KomaMRICore/src/simulation/SimMethods/BlochMagnus/gpu/MagnusQuadKernel.jl

@@ -5,9 +5,9 @@
     M_xy::AbstractVector{Complex{T}}, M_z,
     @Const(p_x), @Const(p_y), @Const(p_z), @Const(p_ΔBz), @Const(p_T1), @Const(p_T2), @Const(p_ρ), N_spins,
     @Const(s_Gx), @Const(s_Gy), @Const(s_Gz), @Const(s_Δt), @Const(s_Δf), @Const(s_B1), @Const(s_ψ), @Const(s_ADC), s_length,
-    ::Val{MOTION}, ::Val{USE_WARP_REDUCTION}, ::Val{HAS_ADC},
+    ::Val{MOTION}, ::Val{USE_WARP_REDUCTION}, ::Val{HAS_ADC}, ::Val{SKIP_RELAXATION},
     sim_method::SM
-) where {T, MOTION, USE_WARP_REDUCTION, HAS_ADC, SM<:Union{BlochMagnusQuad2,BlochMagnusQuad4}}
+) where {T, MOTION, USE_WARP_REDUCTION, HAS_ADC, SKIP_RELAXATION, SM<:Union{BlochMagnusQuad2,BlochMagnusQuad4}}
 
     @uniform N = @groupsize()[1]
     i_l = @index(Local, Linear)
@@ -36,9 +36,11 @@
         ΔBz = p_ΔBz[i]
         Mxy_r, Mxy_i = reim(M_xy[i])
         Mz = M_z[i]
-        ρ = p_ρ[i]
-        T1 = p_T1[i]
-        T2 = p_T2[i]
+        if !SKIP_RELAXATION
+            ρ = p_ρ[i]
+            T1 = p_T1[i]
+            T2 = p_T2[i]
+        end
 
         ψ_start = s_ψ[1]
         if !iszero(ψ_start)
@@ -75,14 +77,16 @@
                 sim_method,
             )
             M_norm = mag_norm(T, Mxy_r, Mxy_i, Mz)
-            Mxy_new_r, Mxy_new_i, Mz_new = rotate_magnetization(θx, θy, θz, Mxy_r, Mxy_i, Mz, T)
-            Mxy_new_r, Mxy_new_i, Mz_new = restore_mag_norm(M_norm, Mxy_new_r, Mxy_new_i, Mz_new) # For reduced float precision only.
-
-            E1 = exp(-Δt / T1)
-            E2 = exp(-Δt / T2)
-            Mxy_r = Mxy_new_r * E2
-            Mxy_i = Mxy_new_i * E2
-            Mz = Mz_new * E1 + ρ * (T(1) - E1)
+            Mxy_r, Mxy_i, Mz = rotate_magnetization(θx, θy, θz, Mxy_r, Mxy_i, Mz, T)
+            Mxy_r, Mxy_i, Mz = restore_mag_norm(M_norm, Mxy_r, Mxy_i, Mz) # For reduced float precision only.
+
+            if !SKIP_RELAXATION
+                E1 = exp(-Δt / T1)
+                E2 = exp(-Δt / T2)
+                Mxy_r *= E2
+                Mxy_i *= E2
+                Mz = Mz * E1 + ρ * (T(1) - E1)
+            end
 
             x0, y0, z0 = x1, y1, z1
             Bx_0, By_0, Bz_0 = Bx_1, By_1, Bz_1

KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl

@@ -44,6 +44,8 @@ Base.view(p::PreallocResult, i::UnitRange) = p
 
 """Default preallocation function."""
 prealloc(sim_method::SimulationMethod, backend::KA.Backend, obj::Phantom{T}, M::Mag{T}, max_block_length::Integer, groupsize) where {T<:Real} = DefaultPrealloc{T}()
+prealloc(sim_method::SimulationMethod, backend::KA.Backend, obj::Phantom, M::Mag, max_block_length::Integer, groupsize, ::Val) =
+    prealloc(sim_method, backend, obj, M, max_block_length, groupsize)
 
 include("BlochSimple/BlochSimple.jl")
 include("Bloch/cpu/BlochCPU.jl")