Mentions légales du service

Skip to content

Modify_mf

FANG Chengran requested to merge modify_mf into master

This pull request focuses on improving the performance of the Matrix Formalism module.

Main changes:

1. If all permeabilities are zero, the Matrix Formalism module separately computes the eigenfunctions for each compartment.

  • If all permeabilities are zero, lap_eig becomes a indexed structure. lap_eig(icmpt) contains the eigvalues, eigfuncs, totaltime and length_scales for compartment icmpt. compute_laplace_eig computes the eigendecomposition for each compartment separately if all permeabilities are zero.
- lap_eig = compute_laplace_eig(femesh, setup.pde, eiglim, setup.mf.neig_max);
+ lap_eig = compute_laplace_eig(femesh, setup.pde, setup.mf, savepath);
  • mf_jn becomes a cell of size 1 x ncompartment. mf_jn{icmpt} contains the jn for compartment icmpt.
- mf_jn = compute_mf_jn(lap_eig.values, setup);
+ mf_jn = compute_mf_jn(lap_eig, setup);
  • diffusion_tensor becomes an array of size [3 x 3 x nsequence x ncompartment]. diffusion_tensor_all is the diffusion tensor for the whole geometry. a is a 1 x ncompartment cell containing the first-order moment.
- diffusion_tensor = compute_mf_diffusion_tensor(femesh, lap_eig, mf_jn);
+ [diffusion_tensor, diffusion_tensor_all, mf_jn, a] = compute_mf_diffusion_tensor(femesh, setup, lap_eig);

2. Improve the computational efficiency of the Matrix Formalism module.

  • Replace expm with expmv to avoid the time-and-memory-consuming matrix exponential operation.
  • Change the default sigma parameter of eigs. The new setting can bring a ~5x speed-up without degrading the precision.
- [funcs, values] = eigs(K+Q, M, neig_max, "smallestreal", otherparams{:});
+ [funcs, values] = eigs(K+Q, M, neig_max, -1e-8, otherparams{:});
  • Remove parfor because the time-consuming operations, such as eigs, eig, expm, expmv, and matrix multiplication, are multithreading and can use all available cores. parfor only adds unnecessary overhead.
  • TODO: Support multiple GPU.

3. Implement the surface relaxation approximation. https://github.com/jingrebeccali/SpinDoctor/pull/3/commits/261b247424835c6f5d6af004fff3e07c0c769eb0

  • Grebenkov treats the transport phenomenon across the cell membranes as surface relaxation in his paper Laplacian Eigenfunctions in NMR. I. A Numerical Tool. Based on this idea, there is an alternative way to implement permeable interfaces (see the section Alternative Implementation of Uniform Surface Relaxation in Grebenkov's paper). The advantage of the surface relaxation approximation is that it only requires the eigendecomposition with the zero Neumann boundary condition to conduct simulation with permeable interfaces.
  • The surface relaxation approximation degrades the precision. But its precision is still acceptable (relative error <= 1% using BTPDE as the reference) with appropriate truncation of the eigendecomposition.
  • I add a new option surf_relaxation to the setup.mf.
setup.mf.surf_relaxation = true;

4. Optimize the I/O module for saving and loading simulated results.

5. Add support for PLY and OFF type meshes. https://github.com/jingrebeccali/SpinDoctor/pull/3/commits/3d15cc2eb3d3b1309183aab8c987c7e289a80634

Minor changes:

Merge request reports