Workshop on fast solvers for fractional diffusion problems
UPDATE: 16 March, 2020
Please note that, in light of the covid19 situation, this workshop will not go ahead as a facetoface meeting. updates will be sent to all registered attendees.
2 April, 2020
9:3017:00
Technology and Innovation Centre,
University of Strathclyde, Glasgow, G1 1RD
This workshop will bring together those actively working on solvers for fractional diffusion problems, and researchers whose tools and expertise may help to solve some of the remaining challenges in the field.
Invited speakers:

Lehel Banjai (HeriotWatt) [abstract]

Stefano Massei (EPFL, Eindhoven) [abstract]

Mariarosa Mazza (Insubria) [abstract]

Spyros Pougkakiotis (Edinburgh) [abstract]

Andy Wathen (Oxford) [abstract]
The programme and abstracts are below.
There will also be a conference dinner on the evening of April 2. More details will be available shortly.
There is no registration fee, but for catering purposes please register using the link below.
For additional information, contact Jennifer Pestana
Abstracts
Spectral Fractional Diffusion via Decoupled Reaction Diffusion Equations
Lehel Banjai
Joint work with J. Markus Melenk and Ch. Schwab
We discuss the discretization of the spectral fractional Laplacian and the efficient solution of the arising linear systems. We consider two representations of the fractional Laplacian  one as the DirichlettoNeumann operator of a degenerate but local elliptic problem posed in an infinite cylinder (sometimes called the CaffarelliSilvestre extension) and the second the representation by the DunfordTaylor integral. We show how a combination of an hpFEM in the extended direction and a diagonalization procedure leads to a decoupled system of reaction diffusion systems for the first approach. A careful choice of parameters in the hpFEM is needed in order to obtain stable results. A quadrature of the DunfordTaylor integral leads directly to such a system as well. We describe how this decomposition into reaction diffusion problems aids both the analysis, e.g., of a fully hpFEM discretization, and the efficient solution. Finally we compare the two approaches.
Parallel solution of large sparse systems by direct and hybrid methods
Iain Duff
We discuss algorithms and codes for the solution of sparse systems that we have developed in an EU Horizon 2020 Project called NLAFET.
We used two approaches to get good single node performance. For symmetric sys tems we used taskbased algorithms based on an assembly tree representation of the factorization. We then used runtime systems for scheduling the computation on both multicore CPU nodes and GPU nodes [4]. The second approach was to design a new parallel threshold Markowitz algorithm [2]. This second approach was designed for highly unsymmetric systems, and we do not consider this further in this talk.
For the first approach, in addition to discussing the matrix factorization, we consider the efficient solution of systems using the factors from this factorization [1].
We then extended the scope of these two approaches to exploit distributed memory parallelism. In both cases, we partition the matrix and run our singlenode code on each partition. We discuss the first case, where we base our work on the block Cimmino algorithm [3] using the ABCD software package coded by Zenadi in Toulouse [6]. We partition the matrix using a numerically aware partitioning algorithm [5]. The kernel for this algorithm is the direct factorization of a symmetric indefinite submatrix for which we use the symmetric treebased code.
We show the performance of our codes on industrial strength large test problems on a heterogeneous platform. Our codes that are available on github are shown to outperform other stateoftheart codes.
Rational Krylov for fractional diffusion problems: convergence and pole selection
Stefano Massei
Evaluating the action of a matrix function on a vector, that is x = f(M)v, is an ubiquitous task in applications. When the matrix M is large, subspace projection methods, such as the rational Krylov method, are usually employed. In this work, we provide a quasioptimal pole choice for rational Krylov methods when f is either CauchyStieltjes or LaplaceStieltjes for a positive definite matrix M. This includes the case f(z) = z^(−α) with α > 0, which arises when solving fractional diffusion problems.
Then, we consider the case when the argument M has the Kronecker structure M = I ⊗ A + B ⊗ I, and v is obtained by vectorizing a lowrank matrix. This finds application, for instance, in solving two dimensional FDEs on rectangular domains. We introduce an error analysis for the numerical approximation of x. Pole choices and explicit convergence bounds are given also in this case.
Multigrids for both isotropic and anisotropic spacefractional diffusion equations
Mariarosa Mazza
Joint work with Marco Donatelli, Rolf Krause and Ken Trotti
We are interested in a finite difference discretization of a twodimensional timedependent space fractional diffusion equation (FDE) and in the solution of the related linear systems. Our focus is both on the case where the fractional orders are close to each other (isotropic case), and the one where they are not (anisotropic case). Driven by the fact that the discretization matrices have a blockToeplitzToeplitzblockslike structure, and by the wellknown negative results on multilevel circulant preconditioning, we opt for a multigrid approach. In this framework, the spectral properties of the matrices and the isotropic/anisotropic nature of the problem guide the selection of both projector and smoother.
In the isotropic case, based on the spectral analysis of the coefficient matrices, we define a multigrid method with classical linear interpolation as grid transfer operator and dampedJacobi as smoother. In the anisotropic case, inspired by certain multigrid strategies for integer order anisotropic diffusion equations given in literature, we replace the classical linear interpolation with a semicoarsening technique. Moreover, we estimate the Jacobi relaxation parameter by using an automatic spectral based procedure. A further improvement in the robustness of the proposed method with respect to the anisotropy of the problem is attained employing the resulting Vcycle with semicoarsening as smoother inside an outer fullcoarsening.
Fast Solutions Methods for Convex Quadratic Fractional Differential Equation Optimization
Spyridon Pougkakiotis
Joint work with John W. Pearson, Santolo Leveque and Jacek Gondzio
We present an optimization method suitable for solving convex quadratic Partial Differential Equation (PDE) optimization problems with box constraints on the state and control variables. In particular, we assume that we are given an arbitrary PDE constrained inverse problem, an associated discretization method, and that the resulting sequences of discrete matrices belong to the class of Generalized Locally Toeplitz (GLT) sequences. Then, we propose the use of an Alternating Direction Method of Multipliers (ADMM) for solving the discretized optimization problems. We show that the linear systems required to be solved during the iterations of ADMM preserve the GLT structure of the initial problem matrices. Using this structure, we present and analyze some general methodologies for preconditioning efficiently such linear systems, and solving them using an appropriate Krylov subspace method.
Subsequently, we focus on a certain class of convex quadratic optimization problems with Fractional Differential Equation (FDE) constraints. Discretized versions of FDEs involve large dense linear systems. In order to overcome this difficulty, we design a recursive linear algebra, based on the fast Fourier trans form (FFT). Focusing on timedependent 2dimensional FDEs, we manage to keep the storage requirements linear, with respect to the grid size N, while en suring an order N log N computational complexity per iteration of the Krylov solver inside ADMM. We implement the proposed method, and demonstrate its robustness and efficiency, through a series of experiments over different setups of the FDE optimization problem.
Parallel preconditioning for timedependent PDEs and PDE control
Andy Wathen
Joint work with Elle McDonald, Jennifer Pestana, Anthony Goddard and Federico Danieli
We present a novel approach to the solution of timedependent PDEs via the socalled monolithic or allatonce formulation. This approach will be explained for simple parabolic and hyperbolic problems and its utility in the context of PDE constrained optimization problems will be elucidated.
The underlying linear algebra includes circulant matrix approximations of Toeplitzstructured matrices and allows for effective parallel implementation. Supporting theoretical results and simple computational results will be shown for the heat equation and the wave equation which indicate the potential as a parallelintime method.