Edit on GitHub

Matrix Formalism

In this example we will consider the matrix formalism approach for a geometry of cylinders.

using SpinDoctor
using LinearAlgebra
using Printf

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

setup = DiskSetup(;
    ncell = 5,
    nsidewall = 20,
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.3,
    layersizes = [0.3, 0.6, 1.0],
    ecs = ConvexHullECS(; margin = 2.0),
    refinement = 0.2,
)
DiskSetup{Float64, ConvexHullECS{Float64}}
  ncell: Int64 5
  layersizes: Array{Float64}((3,)) [0.3, 0.6, 1.0]
  rmin: Float64 2.0
  rmax: Float64 6.0
  dmin: Float64 0.2
  dmax: Float64 0.3
  nsidewall: Int64 20
  ecs: ConvexHullECS{Float64}
  refinement: Float64 0.2

We also define coefficients for the three cell compartments and the ECS.

nlayer = length(setup.layersizes)
coeffs = coefficients(
    setup;
    D = (; cell = [0.001 * I(2) for _ = 1:nlayer], ecs = 0.002 * I(2)),
    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.0, nlayer),
        cell_ecs = 1e-5,
        ecs = 0.0,
    ),
    γ = 2.67513e-4,
)
(D = [[0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.001 0.0; 0.0 0.001], [0.002 0.0; 0.0 0.002]], T₂ = [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf], κ = [0.0, 0.0, 0.0, 0.0, 0.0001, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], ρ = ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im], γ = 0.000267513)

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

mesh, surfaces, cells = create_geometry(setup; recreate = true)
plot_mesh(mesh, 1:3)
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);

We may also compute some useful quantities, including a scalar diffusion coefficient from the diffusion tensors.

volumes = get_cmpt_volumes(model.mesh)
D_avg = 1 / 2 * tr.(model.D)' * volumes / sum(volumes)
0.0016011181178396726

The eigenfunctions of the diffusive part of the Bloch-Torrey operator forms a good basis for the finite element function space. The basis may be truncated at a certain level, thus reducing the number of degrees of freedom. We here opt for 400 eigenfunctions.

laplace = Laplace(; model, matrices, neig_max = 200)
lap_eig = solve(laplace)
(values = [1.915839291004178e-17, 6.7118863591245674e-6, 8.438286101062178e-6, 1.0084227858552708e-5, 1.4621351889353833e-5, 1.8277723969996635e-5, 2.7423920424831355e-5, 3.23322804717178e-5, 5.4280133622259764e-5, 7.740820775138158e-5  …  0.004758195281941791, 0.0047675037325057065, 0.0048476024682635325, 0.004862775015549028, 0.004909989265933959, 0.004921768502855099, 0.004933043734726779, 0.0049728782513586214, 0.004973368971373364, 0.0049939415965742475], funcs = [-0.043085376661377565 -0.08021114748991558 … 0.0004725399482766485 -0.0009582518514917525; -0.04308537666137675 -0.08023582916270051 … -1.976032224183864e-5 -0.0003250879082174944; … ; -0.043085376661390964 0.0023655554154393688 … -0.0011897510228990296 -0.049964545336295185; -0.04308537666139181 0.0028903778870446304 … -0.0009332061934443111 -0.07113525264648161], moments = [[0.27408871947630864 0.06674864871057724 … -0.05011076435675932 0.01485253843874218; 0.06674864871057752 -0.6947568495379004 … 0.09589339483302907 0.019286277610647124; … ; -0.05011076435675935 0.09589339483302904 … -2.0409933219680734 -0.01228149558980211; 0.014852538438742097 0.019286277610647145 … -0.012281495589802105 -2.0029963235400126], [0.3196485474022607 5.320098530099002 … -0.02171852167091664 -0.08372249789585302; 5.320098530099002 0.6859490678838989 … 0.043359052672653754 -0.028074438776370177; … ; -0.021718521670916667 0.04335905267265379 … -4.773588286751684 -0.20365271482981448; -0.08372249789585304 -0.028074438776370007 … -0.2036527148298144 4.8310485104247]], massrelax = [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

The resulting eigenvalues may be represented as length scales, describing the wavelength of the eigenfunctions.

length_scales = eig2length.(lap_eig.values, D_avg)
200-element Vector{Float64}:
  2.87198403096981e7
 48.522063437907505
 43.2747368454441
 39.5858735941269
 32.87515385500371
 29.403599689540336
 24.004727962449486
 22.10769120862417
 17.062438719604817
 14.287888503824439
  ⋮
  1.8206062855907188
  1.8055023568792916
  1.802683444679709
  1.7939952437693798
  1.7918471791195147
  1.7897982361877658
  1.782615363938742
  1.782527416852779
  1.778852052141997

We may also further truncate the eigenfunction basis, if we are satisfied skipping features below a threshold length scale of 3 micrometers.

length_scale = 1.0
λ_max = length2eig(length_scale, D_avg)
lap_eig = limit_lengthscale(lap_eig, λ_max)
(values = [1.915839291004178e-17, 6.7118863591245674e-6, 8.438286101062178e-6, 1.0084227858552708e-5, 1.4621351889353833e-5, 1.8277723969996635e-5, 2.7423920424831355e-5, 3.23322804717178e-5, 5.4280133622259764e-5, 7.740820775138158e-5  …  0.004758195281941791, 0.0047675037325057065, 0.0048476024682635325, 0.004862775015549028, 0.004909989265933959, 0.004921768502855099, 0.004933043734726779, 0.0049728782513586214, 0.004973368971373364, 0.0049939415965742475], funcs = [-0.043085376661377565 -0.08021114748991558 … 0.0004725399482766485 -0.0009582518514917525; -0.04308537666137675 -0.08023582916270051 … -1.976032224183864e-5 -0.0003250879082174944; … ; -0.043085376661390964 0.0023655554154393688 … -0.0011897510228990296 -0.049964545336295185; -0.04308537666139181 0.0028903778870446304 … -0.0009332061934443111 -0.07113525264648161], moments = [[0.27408871947630864 0.06674864871057724 … -0.05011076435675932 0.01485253843874218; 0.06674864871057752 -0.6947568495379004 … 0.09589339483302907 0.019286277610647124; … ; -0.05011076435675935 0.09589339483302904 … -2.0409933219680734 -0.01228149558980211; 0.014852538438742097 0.019286277610647145 … -0.012281495589802105 -2.0029963235400126], [0.3196485474022607 5.320098530099002 … -0.02171852167091664 -0.08372249789585302; 5.320098530099002 0.6859490678838989 … 0.043359052672653754 -0.028074438776370177; … ; -0.021718521670916667 0.04335905267265379 … -4.773588286751684 -0.20365271482981448; -0.08372249789585304 -0.028074438776370007 … -0.2036527148298144 4.8310485104247]], massrelax = [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

Each of the resulting eigenfunctions is represented in the same way as the initial magnetization field ρ.

# eiginds = [2, 3, 5, 10, 15, 20, 30, 40, 50, 100, 150, 200]
ncompartment = length(mesh.points)
fig = Figure()
for i = 1:3, j = 1:4
    ieig = 1 + 4(i - 1) + j
    # ieig = eiginds[4(i - 1) + j]
    ϕ_cmpts = split_field(mesh, lap_eig.funcs[:, ieig])
    ax = Axis(fig[i, j]; title = @sprintf("n = %d, ℓ = %.1f", ieig, length_scales[ieig]))
    ax.aspect = DataAspect()
    for icmpt ∈ 1:ncompartment
        f = mesh.elements[icmpt]
        p = mesh.points[icmpt]
        mesh!(ax, p', f'; color = ϕ_cmpts[icmpt], shading = false)
    end
end
fig

We observe that the first functions have large features, while the higher-index functions have more rapidly varying features. We may now choose a gradient and compute the projection of magnetization field onto the truncated basis.

dir = [1.0, 0.0]
profile = CosOGSE(5000.0, 5000.0, 2)
b = 1000
g = √(b / int_F²(profile)) / coeffs.γ
gradient = ScalarGradient(dir, profile, g)
ScalarGradient{Float64, CosOGSE{Float64}}([1.0, 0.0], CosOGSE{Float64}(5000.0, 5000.0, 2), 4.201554156121445)

The matrix formalism problem is solved in the same way as the BTPDE. The time profile is approximated on 500 points, since it is non-constant.

mf = MatrixFormalism(; model, matrices, lap_eig)
ξ = solve(mf, gradient; ninterval = 500)
2443-element Vector{ComplexF64}:
  0.5563477589496536 + 0.10471282163559548im
  0.5554755231992877 + 0.07550262210833913im
  0.5557633334850008 + 0.01788864769898848im
  0.5589389139213448 - 0.044559322713076395im
  0.5613765988643445 - 0.07189353867828172im
  0.5588618659962727 - 0.044404282354914164im
  0.5554694089630381 + 0.014844559144502471im
  0.5552829113752633 + 0.0784020448015998im
  0.5555853114840529 + 0.05277908163790896im
  0.5601629746335248 + 0.027045735972723195im
                     ⋮
   0.364571986988704 - 0.03909378868348569im
  0.3967393631663856 - 0.05533835768511569im
 0.40210808700490164 - 0.05690220121344265im
  0.3796876179547504 - 0.04776438572493397im
   0.367008793346215 - 0.041287733994012955im
  0.3460773897826393 + 0.04567144110311406im
 0.38420420246041037 + 0.06282044004886186im
 0.33474636569304944 + 0.06596535082632055im
  0.3507587062475337 + 0.06706948467992023im

The resulting magnetization field may be plotted.

plot_field(model.mesh, ξ)

This page was generated using Literate.jl.