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(;
    name = "gaze-into-the-orb",
    ncell = 1,
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.3,
    include_in = false,
    in_ratio = 0.6,
    ecs_shape = :no_ecs,
    ecs_ratio = 0.5,
)
coeffs = coefficients(
    setup;
    D = (; in = 0.002 * I(3), out = 0.002 * I(3), ecs = 0.002 * I(3)),
    T₂ = (; in = Inf, out = Inf, ecs = Inf),
    ρ = (; in = 1.0, out = 1.0, ecs = 1.0),
    κ = (; in_out = 1e-4, out_ecs = 1e-4, in = 0.0, out = 0.0, ecs = 0.0),
    γ = 2.67513e-4,
)

(D = StaticArrays.SMatrix{3, 3, Float64, 9}[[0.002 0.0 0.0; 0.0 0.002 0.0; 0.0 0.0 0.002]], T₂ = [Inf], κ = [0.0], ρ = ComplexF64[1.0 + 0.0im], γ = 0.000267513)

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

mesh, = create_geometry(setup; recreate = true)
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, 203, 206, 1, 2  …  193, 195, 198, 201, 204, 207, 208, 209, 210, 211], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  211, 211, 211, 211, 211, 211, 211, 211, 211, 211], [0.37355660869356755, 0.03689607615923174, 0.08246197300048284, 0.08127731636444514, 0.03922317586278093, 0.03928460977433125, 0.14988222818755206, 0.1313095336915274, 0.03689607615923174, 0.23520457062622052  …  0.19376227357860132, 0.12565850044176596, 0.12253756655724848, 0.28303563889884414, 0.3992799893936035, 0.5075508920811607, 0.3454588957490134, 0.525457096457861, 0.3815625974586602, 4.21925820196383], 211, 211), S = sparse([1, 2, 3, 4, 6, 9, 203, 206, 1, 2  …  193, 195, 198, 201, 204, 207, 208, 209, 210, 211], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  211, 211, 211, 211, 211, 211, 211, 211, 211, 211], [0.01216566858747046, -0.001852126615258414, -0.003968022069524017, -0.004679934079216731, -0.0014553461526621244, -0.00010486639257900713, -0.0002952846867315366, 0.0001899114085013689, -0.001852126615258414, 0.007030387575141143  …  -0.003436479689751879, -0.00012944253545744348, 0.000309898115825891, -0.00175045348399633, -0.0034119911769034707, -0.0028859968250207384, -0.001965415195043322, -0.002440846898950219, -0.004166910625490953, 0.02611193050139627], 211, 211), R = sparse(Int64[], Int64[], Float64[], 211, 211), Mx = SparseArrays.SparseMatrixCSC{Float64, Int64}[sparse([1, 2, 3, 4, 6, 9, 203, 206, 1, 2  …  193, 195, 198, 201, 204, 207, 208, 209, 210, 211], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  211, 211, 211, 211, 211, 211, 211, 211, 211, 211], [0.1079122149947972, -0.007703529643281899, 0.009117197697389452, 0.02540066020198205, 0.031795014097018445, 0.037167175135264484, 0.06943778864722851, -0.024007750291252267, -0.007703529643281899, -0.1389596421111435  …  -0.06334307718016158, 0.09952262767613772, 0.03035912920687382, 0.0547857570293611, 0.19520068189934348, -0.379536686186334, 0.46640018070041245, -0.2627370337587083, 0.30884455338280314, 1.1145685402125936], 211, 211), sparse([1, 2, 3, 4, 6, 9, 203, 206, 1, 2  …  193, 195, 198, 201, 204, 207, 208, 209, 210, 211], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  211, 211, 211, 211, 211, 211, 211, 211, 211, 211], [0.052935191905280336, -0.004945919270948211, 0.036669675021424424, -0.0137799733167113, 0.016081439609033535, -0.00015851119368265765, 0.05391267967654271, 0.018090993284902147, -0.004945919270948211, -0.07247453888303566  …  -0.19503281297560981, -0.1188879630914489, -0.12896288035932965, -0.18433803929738707, -0.1881446999561812, -0.40848684670254903, -0.3168342630029501, -0.8848996228144935, -0.6546160750365931, -6.609992430960907], 211, 211), sparse([1, 2, 3, 4, 6, 9, 203, 206, 1, 2  …  193, 195, 198, 201, 204, 207, 208, 209, 210, 211], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  211, 211, 211, 211, 211, 211, 211, 211, 211, 211], [1.2404135245989876, 0.13051354711001012, 0.2712921620651764, 0.2661145604212933, 0.13566501881961474, 0.13466004162541043, 0.42684034217981787, 0.37236372567645315, 0.13051354711001012, 0.8220308781615863  …  -0.4792548081061305, -0.3447422506701349, -0.3432643975543626, -0.17540276330421314, -0.7207002789418244, -0.7099275293705432, -0.369631046925237, -0.25980720691042747, -0.15759061660833623, -6.16375351486863], 211, 211)], Q = sparse(Int64[], Int64[], Float64[], 211, 211), M_cmpts = Any[sparse([1, 2, 3, 4, 6, 9, 203, 206, 1, 2  …  193, 195, 198, 201, 204, 207, 208, 209, 210, 211], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  211, 211, 211, 211, 211, 211, 211, 211, 211, 211], [0.37355660869356755, 0.03689607615923174, 0.08246197300048284, 0.08127731636444514, 0.03922317586278093, 0.03928460977433125, 0.14988222818755206, 0.1313095336915274, 0.03689607615923174, 0.23520457062622052  …  0.19376227357860132, 0.12565850044176596, 0.12253756655724848, 0.28303563889884414, 0.3992799893936035, 0.5075508920811607, 0.3454588957490134, 0.525457096457861, 0.3815625974586602, 4.21925820196383], 211, 211)], S_cmpts = Any[sparse([1, 2, 3, 4, 6, 9, 203, 206, 1, 2  …  193, 195, 198, 201, 204, 207, 208, 209, 210, 211], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  211, 211, 211, 211, 211, 211, 211, 211, 211, 211], [0.01216566858747046, -0.001852126615258414, -0.003968022069524017, -0.004679934079216731, -0.0014553461526621244, -0.00010486639257900713, -0.0002952846867315366, 0.0001899114085013689, -0.001852126615258414, 0.007030387575141143  …  -0.003436479689751879, -0.00012944253545744348, 0.000309898115825891, -0.00175045348399633, -0.0034119911769034707, -0.0028859968250207384, -0.001965415195043322, -0.002440846898950219, -0.004166910625490953, 0.02611193050139627], 211, 211)], Mx_cmpts = Vector{Any}[[sparse([1, 2, 3, 4, 6, 9, 203, 206, 1, 2  …  193, 195, 198, 201, 204, 207, 208, 209, 210, 211], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  211, 211, 211, 211, 211, 211, 211, 211, 211, 211], [0.1079122149947972, -0.007703529643281899, 0.009117197697389452, 0.02540066020198205, 0.031795014097018445, 0.037167175135264484, 0.06943778864722851, -0.024007750291252267, -0.007703529643281899, -0.1389596421111435  …  -0.06334307718016158, 0.09952262767613772, 0.03035912920687382, 0.0547857570293611, 0.19520068189934348, -0.379536686186334, 0.46640018070041245, -0.2627370337587083, 0.30884455338280314, 1.1145685402125936], 211, 211)], [sparse([1, 2, 3, 4, 6, 9, 203, 206, 1, 2  …  193, 195, 198, 201, 204, 207, 208, 209, 210, 211], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  211, 211, 211, 211, 211, 211, 211, 211, 211, 211], [0.052935191905280336, -0.004945919270948211, 0.036669675021424424, -0.0137799733167113, 0.016081439609033535, -0.00015851119368265765, 0.05391267967654271, 0.018090993284902147, -0.004945919270948211, -0.07247453888303566  …  -0.19503281297560981, -0.1188879630914489, -0.12896288035932965, -0.18433803929738707, -0.1881446999561812, -0.40848684670254903, -0.3168342630029501, -0.8848996228144935, -0.6546160750365931, -6.609992430960907], 211, 211)], [sparse([1, 2, 3, 4, 6, 9, 203, 206, 1, 2  …  193, 195, 198, 201, 204, 207, 208, 209, 210, 211], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  211, 211, 211, 211, 211, 211, 211, 211, 211, 211], [1.2404135245989876, 0.13051354711001012, 0.2712921620651764, 0.2661145604212933, 0.13566501881961474, 0.13466004162541043, 0.42684034217981787, 0.37236372567645315, 0.13051354711001012, 0.8220308781615863  …  -0.4792548081061305, -0.3447422506701349, -0.3432643975543626, -0.17540276330421314, -0.7207002789418244, -0.7099275293705432, -0.369631046925237, -0.25980720691042747, -0.15759061660833623, -6.16375351486863], 211, 211)]], G = Any[[0.10103771509203045 0.00037709022465919756 0.7737757938480939; -0.11269559662756468 -0.08688254944144276 0.9861713658375841; … ; 0.0 0.0 0.0; 0.0 0.0 0.0]], volumes = [260.15298925093447])

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
φ = -π / 6
R = [cos(φ) sin(φ) 0; -sin(φ) cos(φ) 0; 0 0 1]
g⃗(t) = 1.0 * R * [sin(2π * t / TE), sin(20π * t / TE) / 5, cos(2π * t / TE)]
gradient = GeneralGradient(; g⃗, TE)

GeneralGradient{Float64, typeof(Main.g⃗)}(Main.g⃗, 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{Float64}()

Plotter{Float64}(1, 1, Observable{Vector{Float64}} with 0 listeners. Value:
Float64[], Observable{Vector{Float64}} with 0 listeners. Value:
Float64[], Observable{Vector{GeometryBasics.Vec{3, Float32}}} with 0 listeners. Value:
GeometryBasics.Vec{3, Float32}[[0.0, 0.0, 0.0]], Observable{Vector{GeometryBasics.Vec{3, Float32}}} with 0 listeners. Value:
GeometryBasics.Vec{3, Float32}[[0.0, 0.0, 0.0]], Observable{Vector{ComplexF64}} with 0 listeners. Value:
ComplexF64[], Observables.Observable{Vector{Float64}}[Observable{Vector{Float64}} with 0 listeners. Value:
Float64[];;], Observables.Observable{Vector{Float64}}[Observable{Vector{Float64}} with 0 listeners. Value:
Float64[];;], Observable{Vector{Float64}} with 0 listeners. Value:
Float64[], 1.0, 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.