Custom gradients

Here we will consider a custom gradient applied to a spherical geometry.

We start by loading SpinDoctor and a Makie plotting backend.

using SpinDoctor
using LinearAlgebra

if haskey(ENV, "GITHUB_ACTIONS")
    using CairoMakie
else
    using GLMakie
end

setup = SphereSetup(; ncell = 1, layersizes = [0.6, 1.0], rmin = 5.0, rmax = 5.0)
nlayer = length(setup.layersizes)
coeffs = coefficients(
    setup;
    D = (; cell = [0.002 * I(3) for _ = 1:nlayer], ecs = 0.002 * I(3)),
    T₂ = (; cell = fill(Inf, nlayer), ecs = Inf),
    ρ = (; cell = fill(1.0, nlayer), ecs = 1.0),
    κ = (;
        cell_interfaces = fill(1e-4, nlayer - 1),
        cell_boundaries = fill(0, nlayer),
        cell_ecs = 1e-4,
        ecs = 0,
    ),
    γ = 2.67513e-4,
)

(D = [[0.002 0.0 0.0; 0.0 0.002 0.0; 0.0 0.0 0.002], [0.002 0.0 0.0; 0.0 0.002 0.0; 0.0 0.0 0.002]], T₂ = [Inf, Inf], κ = [0.0001, 0.0, 0.0], ρ = ComplexF64[1.0 + 0.0im, 1.0 + 0.0im], γ = 0.000267513)

We then proceed to build the geometry and finite element mesh.

mesh, = create_geometry(setup)
plot_mesh(mesh)

The mesh looks good, so we may then proceed to assemble the biological model and the associated finite element matrices.

model = Model(; mesh, coeffs...)
matrices = assemble_matrices(model)

(M = sparse([1, 2, 3, 4, 6, 9, 73, 1, 2, 3  …  344, 345, 144, 145, 337, 340, 342, 343, 344, 345], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  344, 344, 345, 345, 345, 345, 345, 345, 345, 345], [0.3538648858286587, 0.06831793394256867, 0.0680652700824187, 0.0636492943558472, 0.07783503429254535, 0.07599735315527877, 0.17693244291432936, 0.06831793394256867, 0.4506764076833036, 0.08259293269344262  …  0.5838917813020721, 0.06460680280924767, 0.05505495791417643, 0.13194328595449345, 0.05278343476641638, 0.05331629231583808, 0.059218527812798794, 0.04670406176559619, 0.06460680280924767, 0.309084908892378], 345, 345), S = sparse([1, 2, 3, 4, 6, 9, 73, 1, 2, 3  …  344, 345, 144, 145, 337, 340, 342, 343, 344, 345], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  344, 344, 345, 345, 345, 345, 345, 345, 345, 345], [0.007358089279625656, -0.0016350073918630784, -0.0018278920198483363, -0.0019934985583552143, -0.0013949775089865985, -0.00031163541688142163, -0.00019507838369100806, -0.0016350073918630784, 0.007281762873837704, -7.389775186078719e-5  …  0.01154696549931349, -0.002232345462802864, 0.0016867175778228576, -0.002482153735521645, -0.00021478510281502884, -0.0009047876610191673, -0.0027502954558861975, -0.0011820702133078125, -0.002232345462802864, 0.008079720053529858], 345, 345), R = sparse(Int64[], Int64[], Float64[], 345, 345), Mx = SparseArrays.SparseMatrixCSC{Float64, Int64}[sparse([1, 2, 3, 4, 6, 9, 73, 1, 2, 3  …  344, 345, 144, 145, 337, 340, 342, 343, 344, 345], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  344, 344, 345, 345, 345, 345, 345, 345, 345, 345], [0.1735612784593014, 0.001825118025305791, 0.017060435080313293, 0.03309534053532403, 0.06204314579605384, 0.07285245164489115, 0.07208755390823761, 0.001825118025305791, -0.204875043738325, -0.02304651829112287  …  -0.2846696657687894, 0.0034043660044780415, 0.011412780344337344, 0.06759617651435008, 0.056053634439132506, 0.04564055644917018, 0.03349412329585898, 0.013256193877592004, 0.0034043660044780415, 0.1539200795168264], 345, 345), sparse([1, 2, 3, 4, 6, 9, 73, 1, 2, 3  …  344, 345, 144, 145, 337, 340, 342, 343, 344, 345], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  344, 344, 345, 345, 345, 345, 345, 345, 345, 345], [-0.004712300185233508, -0.012327267175401762, 0.02605684710257862, -0.028303806256626735, 0.02313654097613725, -0.015630764924537396, -0.0023561500926169883, -0.012327267175401762, -0.17894134313625107, 0.016176676611297097  …  0.23781989740653675, 0.015149746787874235, 0.023534664552304036, -0.013906809918190655, 0.004613809743999906, -0.02210501542956199, 0.02633377646333982, -0.021477217581959446, 0.015149746787874235, 0.0033495015039573597], 345, 345), sparse([1, 2, 3, 4, 6, 9, 73, 1, 2, 3  …  344, 345, 144, 145, 337, 340, 342, 343, 344, 345], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  344, 344, 345, 345, 345, 345, 345, 345, 345, 345], [0.8336932850948208, 0.16417071832443098, 0.161260978693768, 0.14773261617940883, 0.1751443181188871, 0.16604367664110492, 0.32960911861048403, 0.16417071832443098, 1.0419841434008632, 0.19411262060021367  …  -2.5388254106763566, -0.29128227647519944, -0.22317478025564175, -0.5485787559714167, -0.23597706777210894, -0.24043059445352288, -0.2645088214053576, -0.2139758290967012, -0.29128227647519944, -1.3933884181498695], 345, 345)], Q = sparse([1, 2, 3, 4, 6, 9, 74, 75, 76, 77  …  69, 70, 71, 72, 137, 140, 142, 143, 144, 145], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  145, 145, 145, 145, 145, 145, 145, 145, 145, 145], [6.13665806570885e-5, 1.187138463239236e-5, 1.1807457521204374e-5, 1.0993365821940555e-5, 1.3483643161544254e-5, 1.321072952000696e-5, -6.13665806570885e-5, -1.187138463239236e-5, -1.1807457521204374e-5, -1.0993365821940555e-5  …  -1.099336582194068e-5, -1.1807457521204401e-5, -1.1871384632392299e-5, -6.136658065708865e-5, 1.3210729520007053e-5, 1.3483643161544223e-5, 1.099336582194068e-5, 1.1807457521204401e-5, 1.1871384632392299e-5, 6.136658065708865e-5], 345, 345), M_cmpts = Any[sparse([1, 2, 3, 4, 6, 9, 73, 1, 2, 3  …  64, 65, 66, 67, 68, 69, 70, 71, 72, 73], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  73, 73, 73, 73, 73, 73, 73, 73, 73, 73], [0.3538648858286587, 0.06831793394256867, 0.0680652700824187, 0.0636492943558472, 0.07783503429254535, 0.07599735315527877, 0.17693244291432936, 0.06831793394256867, 0.4506764076833036, 0.08259293269344262  …  0.22082938354572504, 0.22815866117610273, 0.22557869831566557, 0.22443850481826846, 0.2327945318355359, 0.21194685947549427, 0.22770951417099286, 0.22533820384165035, 0.17693244291432944, 10.399433828588265], 73, 73), sparse([1, 2, 3, 4, 6, 9, 73, 74, 75, 76  …  271, 272, 71, 72, 264, 267, 269, 270, 271, 272], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  271, 271, 272, 272, 272, 272, 272, 272, 272, 272], [0.9972023140506708, 0.09645170767649658, 0.10723698304697732, 0.06749527162630904, 0.0782842227481457, 0.09401536946205338, 0.12191877977621604, 0.12945412531403466, 0.13989829377624763, 0.11177868745566846  …  0.5838917813020721, 0.06460680280924767, 0.05505495791417643, 0.13194328595449345, 0.05278343476641638, 0.05331629231583808, 0.059218527812798794, 0.04670406176559619, 0.06460680280924767, 0.309084908892378], 272, 272)], S_cmpts = Any[sparse([1, 2, 3, 4, 6, 9, 73, 1, 2, 3  …  64, 65, 66, 67, 68, 69, 70, 71, 72, 73], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  73, 73, 73, 73, 73, 73, 73, 73, 73, 73], [0.007358089279625656, -0.0016350073918630784, -0.0018278920198483363, -0.0019934985583552143, -0.0013949775089865985, -0.00031163541688142163, -0.00019507838369100806, -0.0016350073918630784, 0.007281762873837704, -7.389775186078719e-5  …  -0.0003498322644006131, -0.00040055012128052854, -0.0003958567716487311, -0.0003806412779614011, -0.00039420907493098835, -0.00035588267338141906, -0.00038673154223432766, -0.00038557762878865744, -0.00019507838369101066, 0.025021143529916995], 73, 73), sparse([1, 2, 3, 4, 6, 9, 73, 74, 75, 76  …  271, 272, 71, 72, 264, 267, 269, 270, 271, 272], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  271, 271, 272, 272, 272, 272, 272, 272, 272, 272], [0.01345732989680814, -0.003189259912632434, -0.003257843268274944, -0.002041611880021284, -0.0016921551574630759, -0.0007478926436754853, -0.0008574219106319474, 8.6513511361278e-5, 5.3717703757219466e-5, -0.0008862893941823996  …  0.01154696549931349, -0.002232345462802864, 0.0016867175778228576, -0.002482153735521645, -0.00021478510281502884, -0.0009047876610191673, -0.0027502954558861975, -0.0011820702133078125, -0.002232345462802864, 0.008079720053529858], 272, 272)], Mx_cmpts = Vector{Any}[[sparse([1, 2, 3, 4, 6, 9, 73, 1, 2, 3  …  64, 65, 66, 67, 68, 69, 70, 71, 72, 73], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  73, 73, 73, 73, 73, 73, 73, 73, 73, 73], [0.1735612784593014, 0.001825118025305791, 0.017060435080313293, 0.03309534053532403, 0.06204314579605384, 0.07285245164489115, 0.07208755390823761, 0.001825118025305791, -0.204875043738325, -0.02304651829112287  …  0.2287667176403321, -0.24886164031096694, 0.11259680757951285, 0.06259116708826558, -0.18843230501998737, 0.1781165644742895, -0.09188820405534225, -0.011277398098810214, 0.05428593562996631, 0.0034821708172017712], 73, 73), sparse([1, 2, 3, 4, 6, 9, 73, 74, 75, 76  …  271, 272, 71, 72, 264, 267, 269, 270, 271, 272], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  271, 271, 272, 272, 272, 272, 272, 272, 272, 272], [0.5782561480266276, -0.011859209596281683, 0.02040308638125911, 0.053134673023871114, 0.09134963030799863, 0.1314337939081525, 0.07061986993572975, -0.014172540999988908, 0.019761437308166215, 0.061524030369696045  …  -0.2846696657687894, 0.0034043660044780415, 0.011412780344337344, 0.06759617651435008, 0.056053634439132506, 0.04564055644917018, 0.03349412329585898, 0.013256193877592004, 0.0034043660044780415, 0.1539200795168264], 272, 272)], [sparse([1, 2, 3, 4, 6, 9, 73, 1, 2, 3  …  64, 65, 66, 67, 68, 69, 70, 71, 72, 73], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  73, 73, 73, 73, 73, 73, 73, 73, 73, 73], [-0.004712300185233508, -0.012327267175401762, 0.02605684710257862, -0.028303806256626735, 0.02313654097613725, -0.015630764924537396, -0.0023561500926169883, -0.012327267175401762, -0.17894134313625107, 0.016176676611297097  …  -0.12375248083981746, -0.09142341057519789, 0.22062103595488936, -0.22058226841291323, 0.10195822742870021, 0.04434829604160363, -0.09519508457946632, 0.1032224272760291, -0.04748898887648852, -0.0013734628119910636], 73, 73), sparse([1, 2, 3, 4, 6, 9, 73, 74, 75, 76  …  271, 272, 71, 72, 264, 267, 269, 270, 271, 272], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  271, 271, 272, 272, 272, 272, 272, 272, 272, 272], [0.024702429217985628, -0.028428246868137594, 0.06008897322552349, -0.04067848787214139, 0.03448083248569535, -0.020559947334910005, -0.0016874084686180641, -0.024215582562740465, 0.0660958947362815, -0.06029693525970833  …  0.23781989740653675, 0.015149746787874235, 0.023534664552304036, -0.013906809918190655, 0.004613809743999906, -0.02210501542956199, 0.02633377646333982, -0.021477217581959446, 0.015149746787874235, 0.0033495015039573597], 272, 272)], [sparse([1, 2, 3, 4, 6, 9, 73, 1, 2, 3  …  64, 65, 66, 67, 68, 69, 70, 71, 72, 73], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  73, 73, 73, 73, 73, 73, 73, 73, 73, 73], [0.8336932850948208, 0.16417071832443098, 0.161260978693768, 0.14773261617940883, 0.1751443181188871, 0.16604367664110492, 0.32960911861048403, 0.16417071832443098, 1.0419841434008632, 0.19411262060021367  …  -0.3272824549830587, -0.3396974097065934, -0.3478722467806024, -0.3581737655859464, -0.3831726307562445, -0.357604816551345, -0.4082503649740462, -0.4130175156929712, -0.32960911861048453, -9.06219543850284e-15], 73, 73), sparse([1, 2, 3, 4, 6, 9, 73, 74, 75, 76  …  271, 272, 71, 72, 264, 267, 269, 270, 271, 272], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  271, 271, 272, 272, 272, 272, 272, 272, 272, 272], [3.5311326173590465, 0.32610308774392044, 0.3600328028632884, 0.21712259568421552, 0.24228028319868378, 0.2867545596928539, 0.5166430334232387, 0.5341727215725568, 0.5713637199099459, 0.43889085478170864  …  -2.5388254106763566, -0.29128227647519944, -0.22317478025564175, -0.5485787559714167, -0.23597706777210894, -0.24043059445352288, -0.2645088214053576, -0.2139758290967012, -0.29128227647519944, -1.3933884181498695], 272, 272)]], G = Any[[0.25811780109594107 0.0006728397905596947 1.152689191985024; -0.28757502447740285 -0.22099180090099124 1.4598687665917818; … ; 0.1881585648988089 -0.17669636784950732 -1.1526891919850275; 0.0 0.0 0.0], [0.15366093439774375 -0.10221688947872586 0.7292617675422336; -0.0744171407339271 -0.006841490939836499 0.9988421355749963; … ; -0.1660720335457515 0.14783710428110453 -1.5408927591212391; 0.15743659337080107 -0.011724102759403731 -1.2090246778876539]], volumes = [103.99433828588273, 404.1169688448489])

The Bloch-Torrey PDE takes a magnetic field gradient pulse sequence as an input. We may define our custom three dimensional gradient sequence (given in T/m) as a simple Julia function. The echo time TE (given in microseconds) is also needed.

# Magnetic field gradient
TE = 5000.0
φ = -π / 12
R = [cos(φ) sin(φ) 0; -sin(φ) cos(φ) 0; 0 0 1]
gvec(t) = 1.0 * R * [sin(2π * t / TE), sin(20π * t / TE) / 5, cos(2π * t / TE)]
gradient = GeneralGradient(; gvec, TE)

GeneralGradient{Float64, typeof(Main.gvec)}(Main.gvec, 5000.0)

In order to follow the evolution of the solution during time stepping, we add a Plotter to a list of callbacks. Other available callbacks are Printer for showing time stepping information, and VTKWriter for saving the solution time series for visualization in ParaView.

plotter = Plotter()

Plotter{Float64}(1, 1, Observable(0.0), Observable(ComplexF64[]), Scene (800px, 600px):
  0 Plots
  0 Child Scenes)

We may then define the problem and solve for our gradient (with the callback).

btpde = BTPDE(; model, matrices)
ξ = solve(btpde, gradient; callbacks = [plotter])
plotter.fig

MatrixFormalism also accepts GeneralGradients. However, the ADC methods currently only support ScalarGradients.

This page was generated using Literate.jl.