###### Abstract  
We study numerical methods and algorithms for time-dependent fractional-in-space diffusion problems. The considered anomalous diffusion is modelled by the fractional Laplacian \((-\Delta)^{\alpha}\), \(0<\alpha<1\), following the integral definition. Fractional diffusion is non-local, and the finite element method (FEM) discretization in space leads to a dense stiffness matrix. It is well known that numerically solving such non-local boundary value problems is expensive. Difficulties increase significantly when the problem is time-dependent. The aim of the article is to develop computationally efficient methods and algorithms. There are two main features of our approach. Hierarchical semi-separable (HSS) compression is applied for an approximate solution of the arising linear systems. For time discretization, we use an adaptive forward-backward Euler scheme. The properties of the composite algorithm thus obtained are investigated. In particular, the block representation of HSS compression allowed us to upgrade the HSS solver to efficiently handle varying diagonally perturbed transition matrices corresponding to changing time steps. The contribution of the paper is threefold. The methods are completely constructive, which allows for a clearly structured description of the algorithms. A theoretical estimate of the computational complexity is presented. It shows the advantages of the adaptive time stepping in combination with the HSS solver. Theoretical results are supported by representative numerical experiments. Both sequential and parallel scalability and efficiency are analyzed. The presented results provide convincing proof of the concept of the proposed methods and algorithms. f +
Footnote †: journal:  
1
Footnote †: journal:  
1
Footnote 1: Department of Mathematics and Statistics, University of Toronto, Toronto, Ontario, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada, Canada,, Canada, Canada,  
results for different definitions of fractional Laplacian. In [9], a comparison between the spectral and integral definitions can be found. In this work, we will concentrate on the following integral definition of \((-\Delta)^{\alpha}\):  
\[(-\Delta)^{\alpha}u(x)=C(d,\alpha)\ \text{P.V.}\int_{\mathbb{R}^{d}}\frac{u(x)-u(y)} {|x-y|^{d+2\alpha}}dy,\quad\alpha\in(0,1).\]  
Here P.V. stands for principal value, \(d\) is the number of spatial dimensions, \(\alpha\) is the fractional power, and \(C(d,\alpha)\) is the following constant  
\[C(d,\alpha)=\frac{2^{2\alpha}a\Gamma\left(\alpha+\frac{d}{2}\right)}{\pi^{d/2 }\Gamma(1-\alpha)},\]  
where \(\Gamma\) is the gamma function. The fractional Laplacian is a non-local operator, and after applying a suitable discretization in space, the associated boundary value problems are reduced to solving systems of linear algebraic equations with dense matrices. For real-life problems, solving such systems can be computationally very expensive. If a standard LU or Cholesky decomposition is applied, the computational complexity to solve such a system is \(O(N^{3})\), where \(N\) is the number of degrees of freedom in space. For the considered parabolic problem, we need to solve a series of such linear systems, and the complexity becomes \(O(mN^{3})\), where \(m\) is the number of time steps. The aim of this research is to develop a more efficient numerical method with significantly less computational complexity. In this work, we do not impose any special constraints on the geometry of the space computational domain \(\Omega\). Also, in the finite element method (FEM) discretization, the mesh is unstructured by default.