Error during simulation with T2* decay

[Jan-Meyer01]

What happened?

Hello everyone,

while trying to run a simulation with T2* decay (as described in #540) I ran into an strange error that occurred semi-randomly during simulation:

GC error (probable corruption)
Allocations: 110980447 (Pool: 110977746; Big: 2701); GC: 1651
!!! ERROR in jl_ -- ABORTING !!!

I tried different sequences and phantoms, but couldn't really pin it down. The RAM is not overflowing as far as I can tell, however the error seems to occur more often when trying to increase the resolution of the phantoms to simulate e.g. flash/gre sequences.
The following code should reproduce the error (at least it does on my machine):

using KomaMRI

# Define scanner and sequence --> adapt scanner hardware to Prisma 3T
sys = Scanner()
sys.B0 = 3      # 3T
sys.Gmax = 80   # mT/m
sys.Smax = 200  # mT/m

# read example gre sequence and load oversampled phantom
seq = read_seq(joinpath(dirname(pathof(KomaMRI)), "../examples/1.sequences/ge.seq"))
obj = brain_phantom2D(us=(1,2))

# Define simulation parameters
sim_params = KomaMRICore.default_sim_params()
sim_params["gpu"] = false                       # use CPU, because the amount of resources are usually to much for the GPU
sim_params["sim_method"] = BlochWithT2s(100)    # simulate T2star decay with 100 isochromats

# start the simulation and save raw signal as .mrd
raw = simulate(obj, seq, sys; sim_params)
@info """Finished simulation!"""
file = ISMRMRDFile("TestSignal.mrd")
save(file, raw)
@info """Saved simulated signal!"""

The simulation function BlochWithT2s is incorporated as follows in the SimulationMethods.jl file:

# special Bloch method to simulate T2star
struct BlochWithT2s <: SimulationMethod
	Niso::Int
end

export BlochWithT2s

function KomaMRICore.initialize_spins_state(
    obj::Phantom{T}, 
    sim_method::BlochWithT2s # <-------
) where {T<:Real}
    obj = include_T2s(obj, sim_method.Niso) # <-------
    Nspins = length(obj)
    Mxy = zeros(T, Nspins)
    Mz = obj.ρ
    Xt = Mag{T}(Mxy, Mz)
    return Xt, obj
end

function include_T2s(obj, Niso)
	T2prime = (1 ./ obj.T2s .- 1 ./ obj.T2) .^ -1
	prob(Δw) = T2prime ./ (π .* (1 .+ T2prime .^ 2 .* Δw .^ 2 ))
	Δw_lim = 2π * 200 # -200 Hz to 200 Hz
	obj_aux = copy(obj)
    obj_with_T2s = Phantom()
    for  Δw in range(-Δw_lim, Δw_lim, Niso)
		obj_aux.Δw .= Δw
        obj_aux.ρ .*= prob(Δw) # <-- This is what I meant
		obj_with_T2s += obj_aux
    end
	return obj_with_T2s
end

Any help on tracing the error would be appreciated!

Environment

OS x86_64-linux-gnu
Julia 1.11.6
KomaMRIPlots 0.9.2
KomaMRIFiles 0.9.1
KomaMRI 0.9.1
KomaMRICore 0.9.1
KomaMRIBase 0.9.1

[Jan-Meyer01]

Here are some of the actual sequences I used in case it's helpful (though the example above should also work):
ep2d_slice1_FOV200x200_Res64x64_FA5_TE13_TR30.seq.txt
flash_slice1_FOV200x200_Res64x64_FA5_TE13_TR30.seq.txt