Franz Chouly

Jack S Hale

Alexei Lozinski

Stéphane P A Bordas

Olga Barrera

Keywords: Phenomenon Data Model Prediction, Model improvement, Parameter calibration, ... Measurements Modeling Numerical solve Model uncertainties Data uncertainties conceptualization, formal analysis, investigation, methodology, software, validation, visualization, writing -original draft, writing -review & editing conceptualization, formal analysis, investigation, methodology, software, validation, visualization, writing -original draft, writing -review & editing conceptualization, formal analysis, investigation, methodology, software, validation, visualization, writing -original draft, writing -review & editing

. It is of particular importance to balance these errors to avoid wasted computation time and resources.

On a first time, in Chapters 1 to 4, we focus on a family of models based on linear elliptic fractional partial differential equations. In particular, we look at the discretization errors coming from two particular (families of) numerical methods: finite element methods and rational approximation methods.

) quantification.

A posteriori error estimation

In this section we review the literature on a posteriori error estimation techniques for finite element methods. We focus particularly on the implicit/hierarchical a posteriori error estimation methods.

I would like to thank Dr. Nabile Boussaid and Prof. Roland Becker to have accepted to be in my thesis supervision committee. They both played a role in the fact that I am now about to defend a thesis in applied mathematics. Dr. Boussaid, not only you have been a great teacher during my "préparation à l'agrégation" years but you also have introduced the student in search of new horizons I was to whom would become my thesis supervisor. Prof. Becker, for the motivative and fruitful discussions and because you have suggested the study of the Bank-Weiser estimator to one of my thesis supervisor, estimator which is now the central tool of my thesis work.

Then, I thank my thesis supervisors: Prof. Franz Chouly, Prof. Alexei Lozinski, Prof. Stéphane P. A. Bordas, and Dr. Jack S. Hale. I had the chance to be surrounded by very talented researchers and what is most important, great persons.

Franz, you have introduced me to the field of numerical analysis and to error estimation. Your valuable advice on mathematics or on the researcher life, continuous support and dedicated involvement -despite the constantly increasing (Euclidean) distance between us-have been of great help to me.

Alexei, thank you for your warm welcome in the numerical analysis team in Besançon, for all your help with the proofs and for your bright ideas and incredible mathematical intuition.

Stéphane, I thank you for the very nice (and often rambling) discussions about mathematics or life. Thank you also for your very helpful advice on research and for your energy and positivity in all circumstances.

Last but not least, Jack, I would never thank you enough for your help, your patience and almost daily support on all the aspects of my thesis. This experience was new for both of us but you have been a great supervisor from the beginning to the end. Thank you for vii your valuable advice on programming, numerical methods, mechanical engineering and mathematics, your wide range of expertise has been of great help to me. Thank you for the very nice trips, for the beers and good restaurants in Austin.

I would like to thank my colleagues students, PhD students, post-docs and ATER from Besançon, Claire, Lucie, Isabelle, Quentin, Thao, Benjamin, Loris, Marine, Julien and those I forget. It was always a pleasure to meet you each time I came back to Besançon.

A warm thank you to my colleagues from Belval, my office mates Soumi and Hussein and also to Marco (I miss the morning coffee breaks), Paul, Sofia, Antoine, Michal, Arnaud, Saurabh, Chintan, Milad, Lars, Anas. I want to thank particularly the Messins Hugues and Lucie. Hugues thank you for your help, I hope we will continue our collaboration on random maths problems in the future.

Special thanks to the people from the Vortex, thank you very much for providing such a respectful and wholesome place to share our thoughts, discoveries, joys and pains.

Finally, I will switch to French to properly thank my close friends and family. Je voudrais tout d'abord remercier Loupi, Atchoum, Sparrow et Lola pour avoir accompagné mes journées lors du confinement et les autres également.

Je voudrais remercier chaleureusement mes amis de longue date Benjamin et Roxane pour m'avoir accompagné pendant des années et pour être restés proches même en étant éloignés.

Merci à Laura et Nassim pour les parties de jeux de société endiablées et pour le catsitting. Merci à Chantal et Gérard pour leur accueil toujours chaleureux. Merci à mes grand-parents d'avoir tenté de suivre mes aventures (avec plus ou moins de succès) et de m'avoir toujours soutenu.

Introduction

All scientific fields share the same underlying process of acquiring knowledge. This process is based on a feedback loop, sketched in fig. 1. Information on a phenomenon is collected from data gathering. From this information a model, i.e. a system of mathematical equations depending on diverse parameters, is established. Finally, a prediction is drawn from the model and this prediction can be used to improve the model.

However, researchers know that this is a never-ending loop: each step of the above process induces uncertainties. For example, the collection of data is always corrupted with measurements errors, there is no complex enough (manageable) model to capture reality exactly and finally, numerical experiments demand discrete (which means truncated) inputs to perform and therefore, return an approximation of the solution.

Uncertainty can be divided into two categories: aleatory uncertainty and epistemic uncertainty. Epistemic uncertainty arises from a lack of knowledge and can be reduced by collecting more information on the system we are studying. In contrast, aleatory uncertainty is considered as irreducible and inherent to the system. This distinction is not absolute, the choice to classify uncertainty as aleatory or epistemic is a modeling issue [START_REF] Kiureghian | Aleatory or epistemic? Does it matter?[END_REF].

Thus, the result of any scientific study always comes with uncertainties. This is fine as long as we are aware of aleatory uncertainties and are able to quantify epistemic uncertainties. This is what uncertainty quantification is all about: study the impact of these uncertainties on the result of this process. Our interest here is on epistemic uncertainties arising from numerical implementation, "Numerical uncertainties" in fig. 1. In what follows Figure 1: Different sources of uncertainties in the process of numerical simulation.

xiii

INTRODUCTION

we refer to these uncertainties as errors as it is usually the case in the finite element literature.

There are multiple sources of numerical errors, we can cite the discretization errors induced by the use of a numerical scheme (such as e.g. finite element methods) [START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF][START_REF] Verfürth | A Posteriori Error Estimation Techniques for Finite Element Methods[END_REF], the algebraic errors coming from the use of (inexact) iterative solvers (such as e.g. conjugate

Introduction

Error quantification in finite element methods is based on two complementary concepts: a priori and a posteriori error estimation. A priori error estimation gives global upper bounds on the numerical errors and are often used to prove the convergence of uniform refinement methods (in contrast to adaptive refinement methods we will consider here below).

A posteriori error estimation is meant to provide computable approximations to the numerical errors and analysts should be able to use a posteriori error estimators to assess the quality of the simulation they study. Thus, a special care is taken to avoid as much as possible the dependence of a posteriori error estimators on unknown quantities.

The names a priori and a posteriori are transparent and come from the fact that, a priori error estimation, as it might depend on unknown quantities, can be carried out before the numerical solve while a posteriori error estimation takes advantage of the approximated solution and therefore must be performed after the numerical solve.

A posteriori error estimation is a well established tool in finite element discretization methods. Its origin and use in mesh adaptation techniques for finite element methods goes back to the works of Babuška and Rheinboldt [START_REF] Babuška | On the reliability and optimality of the finite element method[END_REF][START_REF] Babuška | Analysis of Optimal Finite-Element Meshes in R 1[END_REF][START_REF] Babuška | A Posteriori Error Analysis of Finite Element Solutions for One-Dimensional Problems[END_REF][START_REF] Babuška | A-posteriori error estimates for the finite element method[END_REF][START_REF] Babuška | Error Estimates for Adaptive Finite Element Computations[END_REF]. A posteriori error estimators are obtained from the derivation of quantities at the level of the domain, of patches of mesh cells or of cells themselves, quantifying the global and/or local discretization error. xv To fix ideas, let us consider a toy problem and an a posteriori error estimator based on local quantities at the level of the cells. Let Ω be an open connected domain with polygonal boundary Γ, H be an Hilbert space of functions defined on Ω and let F : H → R be a linear functional. We consider the following general (elliptic) partial differential equation: find u in

H such that a (u, v) = F (v) ∀v ∈ H, (1)
where a (•, •) : H ×H → R is a bilinear form. Let T h be a mesh on Ω, where h is a parameter encoding the coarseness of the mesh and composed of cells {T }. Let V h H be a finite element space associated to T h . The finite element problem associated to eq. ( 1) is: find

u h in V h such that a (u h , v h ) = F (v h ) ∀v h ∈ V h . (2)
The finite element discretization error is given in a certain norm by u -u h .

A simple example of boundary value problem is the Poisson equation with homogeneous zero Dirichlet boundary conditions on Ω. In this case, the "natural" space in which the solution u is sought is H := H 1 0 (Ω) the Sobolev space of order 1 of functions vanishing on Γ. For (u, v) ∈ H 2 , the bilinear form a is given by a(u, v) := Ω ∇u • ∇v and the linear form F by F (v) := Ω f v, where e.g. f ∈ L 2 (Ω) is the data of the problem. Equation [START_REF] Abd | Neutron radiography study of water absorption in porous building materials: anomalous diffusion analysis[END_REF] becomes: find u in H such that

Ω ∇u • ∇v = Ω f v ∀v ∈ H. (3)
In this particular case, a mesh T h on Ω consists e.g. in a partition of Ω by a set of cells {T }, the parameter h being a measure of the coarseness of the mesh Ω. A finite element space V h H is e.g. the space of continuous, piecewise linear functions on each cell T of the mesh, vanishing on 0 and 1. The finite element discretization of eq. ( 3) is then: find u h in V h such that

Ω ∇u h • ∇v h = Ω f v h ∀v h ∈ V h . (4)
Thus, the finite element discretization error is given, for example, in the norm of the energy by u -u h := a(u -u h , u -u h ) = ∇u -∇u h L 2 .

A typical a priori error estimation result would, for example, give a bound like the following u -

u h 1 Ch δ u 2 , (5)
where • 1 and • 2 are two different norms, C is an eventual unknown constant and δ is a real number giving the rate of convergence of the method as h → 0 (see e.g. [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]Chapter 3.2]). Notice that in eq. ( 5), the bound depends on the unknown solution u.

In contrast, an a posteriori error estimation depends on known quantities only: the problem data F , the bilinear form a and the finite element discretization u h , and is therefore fully computable. Let η be an a posteriori error estimator, derived from local contributions, such that η 2 :=

T ∈T h η 2 T . (6)
Among the desirable properties for η, we can cite:

xvi INTRODUCTION i) global reliability, the existence of a constant C > 0 independent from the solution u, from u h and from the fineness of T h , such that u -u h Cη + osc, [START_REF] Acosta | A Fractional Laplace Equation: Regularity of Solutions and Finite Element Approximations[END_REF] where osc are the data oscillations, measuring the discretization error on the data of the problem (indeed, the linear form F might depend on quantities that must be discretized in order to solve eq. ( 2) numerically e.g. the data f ∈ L 2 (Ω) in eq. ( 4)) and are usually higher order terms. Often, the constant C is difficult to evaluate. When C is known (equivalently, when C = 1), the a posteriori error estimator is called guaranteed estimator.

ii) local efficiency, the existence of a constant c > 0 independent from u, u h and from the fineness of T h , such that

cη T u |ω T -u h|ω T ω T + osc |ω T , (8)
where ω T is often a patch of cells containing T and osc |ω T is a localization of the data oscillations to this patch. The constant c is usually unknown.

iii) asymptotic exactness, the estimator tends to the exact error as we refine the mesh

η u -u h -→ 1, (9)
the ratio η/ u -u h is called efficiency index.

iv) robustness, the constants C, c and the efficiency index do not depend on the parameters of the problem (e.g. coefficients in a or the finite element space V h ), v) low computational cost, the computation of η should be significantly less expensive than the computation of the finite element solution u h .

Optimally, an estimator should satisfy eqs. [START_REF] Acosta | A Fractional Laplace Equation: Regularity of Solutions and Finite Element Approximations[END_REF] and [START_REF] Acosta | A short FE implementation for a 2d homogeneous Dirichlet problem of a fractional Laplacian[END_REF] with constants C and c as close as possible to 1. This guarantees two things, first the estimator is a good approximation to the exact error η u -u h , [START_REF] Ainsworth | The performance of Bank-Weiser's error estimator for quadrilateral finite elements[END_REF] and second, the estimator captures well the distribution of the error across the mesh T h , which is an essential feature for adaptive mesh refinement algorithms. When not explicitly specified, we consider in the following that the norm • in Properties i) to iii) is the energy norm of the associated problem. The literature on a posteriori error estimation is vast and to give an exhaustive list of a posteriori error estimators would be a colossal task. Our goal here is to give an overview of the principal families of estimators (keeping in mind that these families do not constitute a partition of the set of estimators) and some of the seminal papers associated to their introduction and to the study of Properties i) to v), mostly in the context of linear elliptic PDEs. We refer the reader to the following works for a broader view on the topic [START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF][START_REF] Carstensen | Estimator Competition for poisson Problems[END_REF][START_REF] Carstensen | Axioms of adaptivity[END_REF][START_REF] Ern | Finite Elements I[END_REF][START_REF] Ern | Finite Elements II[END_REF][START_REF] Ern | Finite Elements III[END_REF][START_REF] Neittaanmäki | Reliable Methods for Computer Simulation: Error Control and Posteriori Estimates[END_REF][START_REF] Repin | A Posteriori Estimates for Partial Differential Equations[END_REF][START_REF] Verfürth | A review of a posteriori error estimation and adaptive mesh-refinement techniques[END_REF][START_REF] Verfürth | A Posteriori Error Estimation Techniques for Finite Element Methods[END_REF].

xvii

Explicit residual estimators

Residual error estimators are based on the residuals associated to the discretization of a boundary value problem. The residuals are quantities measuring the misfit of the finite element solution in the weak formulation eq. ( 1) and are derived from the following quantity

R(v) := a (u, v) -a (u h , v) = F (v) -a (u h , v) . (11)
By definition the residual R is a global quantity so a first step is often to localize it to the level of patches of cells or cells themselves.

Explicit residual estimators are based on norms of the quantity R. They are the first type of estimators to have been introduced, in [START_REF] Babuška | A Posteriori Error Analysis of Finite Element Solutions for One-Dimensional Problems[END_REF][START_REF] Babuška | A-posteriori error estimates for the finite element method[END_REF][START_REF] Babuška | Error Estimates for Adaptive Finite Element Computations[END_REF]. They are among the cheapest estimators, therefore satisfy Property v), but they can sometimes largely overestimate the discretization error [START_REF] Carstensen | Estimator Competition for poisson Problems[END_REF]. Property i) has been established in [START_REF] Babuška | Analysis of the Efficiency of an a Posteriori Error Estimator for Linear Triangular Finite Elements[END_REF] and Property ii) in [START_REF] Verfürth | A posteriori error estimators for the Stokes equations[END_REF] for an explicit residual estimator adapted to Stokes equations. A similar proof in the case of Poisson equations can be found in [START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF]. Explicit residual estimators satisfying Property iv) in the context of reaction-diffusion and convection-diffusion have been introduced in [START_REF] Verfürth | Robust a posteriori error estimators for a singularly perturbed reactiondiffusion equation[END_REF][START_REF] Verfürth | Robust A Posteriori Error Estimates for Stationary Convection-Diffusion Equations[END_REF]. Guaranteed explicit residual estimators satisfying Property i) with C = 1 can be found e.g. in [START_REF] Carstensen | Fully Reliable Localized Error Control in the FEM[END_REF].

Flux recovery estimators

Fluxes of the solution u are usually the weighted gradient of u and are sometimes the primary concern (e.g. for elasticity problems, practitioners are often interested in stresses and strains, quantities derived from the fluxes of u). However, the normal component of the fluxes of the approximation u h is generally discontinuous across the cell boundaries. This discontinuity is not present in the fluxes of the solution u, which induces an approximation error.

Based on reconstructions of fluxes from u h by averaging or smoothing techniques, they are used to give an estimation of the error on the fluxes of the solution u. These estimators are obtained by computing the norm of the discrepancy between the fluxes of the finite element solution u h and the reconstructed ones. Flux reconstruction estimators have been introduced for the first time in [START_REF] Babuška | A Posteriori Error Analysis of Finite Element Solutions for One-Dimensional Problems[END_REF] for one-dimensional problems and have been generalized in [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery (SPR) and adaptive finite element refinement[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity[END_REF]. In [START_REF] Ainsworth | Analysis of the Zienkiewicz-Zhu a-posteriori error estimator in the finite element method[END_REF] the authors have shown that flux reconstruction estimators satisfying the superconvergence of the recovered flux [START_REF]Some superconvergence results in the finite element method[END_REF] satisfy Property iii). In [START_REF] Rodríguez | Some remarks on Zienkiewicz-Zhu estimator[END_REF], Property i) is shown for the flux reconstruction estimator introduced in [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF]. In addition, it is shown that this estimator does not satisfy Property iii) in general. Property i) is established for all the flux recovery methods in [START_REF] Bartels | Each averaging technique yields reliable a posteriori error control in FEM on unstructured grids. Part II: Higher order FEM[END_REF][START_REF] Carstensen | Each averaging technique yields reliable a posteriori error control in FEM on unstructured grids. Part I: Low order conforming, nonconforming, and mixed FEM[END_REF] and Property ii) is proved for all the methods in [START_REF] Carstensen | All first-order averaging techniques for a posteriori finite element error control on unstructured grids are efficient and reliable[END_REF]. Some of the flux recovery estimators are particularly easy to implement at a low computational cost (they satisfy Property v)), see [START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF]Chapter 4].

Implicit residual estimators/Hierarchical estimators

These estimators are derived from the solutions of local auxiliary boundary value problems.

Since the estimator at the core of this thesis belongs to this family, we give more details in INTRODUCTION the next section.

Functional estimators

Derived from duality principles reformulating the boundary value problem into a minimization one, functional estimators have been studied, originally for nonlinear problems, in [START_REF] Repin | A posteriori error estimation for nonlinear variational problems by duality theory[END_REF][START_REF] Repin | A posteriori error estimates for approximate solutions of variational problems with functionals of power growth[END_REF]. Such estimators are not local and consequently can not satisfy Property ii) or be used in adaptive mesh refinement methods. Moreover, these estimators do not satisfy Property iv). However, they satisfy Property i) with C = 1 by construction. These properties are detailed in [START_REF] Repin | A posteriori error estimates for approximate solutions to variational problems with strongly convex functionals[END_REF][START_REF] Repin | A posteriori error estimation for variational problems with uniformly convex functionals[END_REF].

Equilibration estimators

They are based on a reconstruction of discrete fluxes that satisfy a certain equilibration equation. The particularity of these estimators is that they are designed to be a guaranteed upper bound of the error i.e. to satisfy Property i) with C = 1. One of the first equilibration technique used in a posteriori error estimation comes from the work of [START_REF] Ladeveze | Error Estimate Procedure in the Finite Element Method and Applications[END_REF]. An estimator based on local solves in patches around the vertices of the mesh is derived in [START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF]. In [START_REF] Luce | A Local A Posteriori Error Estimator Based on Equilibrated Fluxes[END_REF] the authors have derived an estimator computed via dual mesh techniques that satisfies Property i), Property ii) and Property iii) in some particular sense. Robust estimators, satisfying Property iv) have been introduced in [START_REF] Ern | Polynomial-Degree-Robust A Posteriori Estimates in a Unified Setting for Conforming, Nonconforming, Discontinuous Galerkin, and Mixed Discretizations[END_REF][START_REF] Vohralík | Guaranteed and Fully Robust a posteriori Error Estimates for Conforming Discretizations of Diffusion Problems with Discontinuous Coefficients[END_REF].

Implicit residual and hierarchical error estimators

In contrast to explicit residual a posteriori error estimators, implicit residual estimators are not directly based on certain norms of quantities coming from the residuals but are derived from the solutions to auxiliary problems for which the right-hand side data is based on R from eq. [START_REF] Ainsworth | The influence and selection of subspaces for a posteriori error estimators[END_REF].

Hierarchical estimators are closely related to implicit residual estimators. Like implicit residual estimators, hierarchical estimators are based on auxiliary problems. However, for hierarchical estimators these auxiliary problems are solved on particular discrete spaces defined from two (or more) nested (discrete) spaces, as we will see in this subsection.

Generalities about implicit and hierarchical estimators

The idea at the core of implicit error estimation techniques comes from the fact that the discrepancy u -u h is solution to the following boundary value problem

a (u -u h , v) = R(v) ∀v ∈ H. ( 12
)
For example, in the case of eqs. ( 3) and (4), using an integration by parts formula on each cell of the mesh, eq. [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF] becomes

Ω ∇(u -u h ) • ∇v = T ∈T h T r T + E∈∂T E J E , (13)
where r T and J E are certain functions depending on f and u h we do not explicit here for the sake of simplicity. Then, implicit residual estimators are based on solutions to finite element discretizations of eq. [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF], eventually localized to patches of cells or to cells levels. For example, such a problem at the level of a subdomain ω (ω can be e.g. a patch of cells, a cell...) would read: seek e ω in W ω such that

a ω (e ω , v ω ) = R ω (v ω ) ∀v ω ∈ W ω , (14)
where a ω and R ω are localized (and sometimes modified) versions of the bilinear form a and of the residual R and W ω is a local finite element space. This approach rises several questions such as: in which spaces {W ω } ω do we seek the solutions to these local problems ? Do we set these problems at the cell level or at the level of patches of cells ? What boundary conditions do we choose ? There is no definitive answers to these questions, different error estimators (with different properties) arise depending on the decisions we take to tackle these questions.

For example, when the finite element space W ω is based on a hierarchy of spaces W ω W ω , the resulting estimator is called a hierarchical estimator. Hierarchical estimator are an important class of estimators in the family of implicit residual estimators.

The first implicit residual estimators are due to Babuška and Rheinboldt [START_REF] Babuška | A-posteriori error estimates for the finite element method[END_REF][START_REF] Babuška | On the reliability and optimality of the finite element method[END_REF]. They introduced the solves of local Dirichlet boundary value problems on patches of cells around each vertex of the mesh. Implicit residual estimators have been popularized by the work of Bank and Weiser [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF] where three novel estimators are introduced, based on solves of local Neumann problems at the level of cells themselves. One of the estimators introduced in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF] is at the core of this thesis. In [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF] is shown that the introduced estimators satisfy Property i) and a global version of Property ii), under a binding assumption called saturation assumption. The saturation assumption is crucial in the analysis of hierarchical error estimators, especially to prove Property i).

There exist several definitions of the saturation assumption but a common one is the following: Let V V be two nested finite element spaces (e.g. V is defined on a finer mesh or contains piecewise polynomial functions of higher degree than V ) and u V and u V be the corresponding finite element solutions to eq. ( 2). Then, there exists a constant

0 < β < 1 such that u -u V β u -u V . (15)
In other words, the saturation assumption claims that the use of a finer finite element space must improve the accuracy of the approximation. This assumption is very difficult to assert in practice (since it requires the knowledge of the exact solution u) and it has been proved not to be true even in very simple situations, especially in the preasymptotic regime [START_REF] Dörfler | Small data oscillation implies the saturation assumption[END_REF]. This is why, after [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF], many works tried to remove the saturation assumption from the analysis of hierarchical estimators, see e.g. [START_REF] Achchab | Some remarks about the hierarchicala posteriori error estimate[END_REF][START_REF] Araya | An a posteriori Error Estimate of Hierarchical Type[END_REF][START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF][START_REF] Nochetto | Removing the saturation assumption in a posteriori error analysis[END_REF][START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF]. It is in this vein that our novel study in Chapter 1 proposes a proof of Property i) for the third estimator introduced in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF], without requiring the saturation assumption. In addition, a novel implementation in the finite element software FEniCS is proposed in Chapter 2 and numerical evidences showing that this estimator also satisfy Property v) are provided.

INTRODUCTION

Implicit residual estimators satisfying Property i) and Property ii) without the use of the saturation assumption are proposed in [START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF]. A particular version of the third estimator introduced in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF] is adapted to singularly perturbed reaction-diffusion equations and convection-diffusion equations in [START_REF] Verfürth | Robust a posteriori error estimators for a singularly perturbed reactiondiffusion equation[END_REF] and [START_REF] Verfürth | A posteriori error estimators for convection-diffusion equations[END_REF] respectively. In these two papers, the resulting estimators are proved to satisfy Property i) and ii) and Property iv). This same estimator has been also adapted to linear elasticity equations in [START_REF] Verfürth | A review of a posteriori error estimation techniques for elasticity problems[END_REF] and is proved to satisfy Property i) and ii).

In the case of linear elasticity problems, [START_REF] Araya | An a posteriori Error Estimate of Hierarchical Type[END_REF] proposes a proof of a variant of the saturation assumption and a hierarchical error estimator satisfying Property i) and Property ii). In addition, their estimator satisfy Property iv) with respect to the material parameters.

As it has already been discussed in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF], the choice of boundary conditions for the local problems is one of the main parameters influencing the performance of implicit residual estimators. The works [START_REF] Ainsworth | A posteriori error estimators for second order elliptic systems: Part 1. Theoretical foundations and a posteriori error analysis[END_REF][START_REF] Ainsworth | A posteriori error estimators for second order elliptic systems: Part 2. An optimal order process for calculating self-equilibrating fluxes[END_REF] develop a way to choose the local problems boundary conditions following self-equilibration techniques introduced in [START_REF] Kelly | The self-equilibration of residuals and complementarya posteriori error estimates in the finite element method[END_REF]. The resulting estimator satisfy a global version of Property ii) and Property i) under some variant of the saturation assumption. The question of finding appropriate boundary conditions for the local problems also led to novel estimators based on solves in patches of cells in [START_REF] Carstensen | Fully Reliable Localized Error Control in the FEM[END_REF][START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF][START_REF] Morin | Local problems on stars: A posteriori error estimators, convergence, and performance[END_REF]. In [START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF] the boundary conditions inside the patches of cells are derived from equilibration. In [START_REF] Carstensen | Fully Reliable Localized Error Control in the FEM[END_REF] a modification of the estimator derived in [START_REF] Babuška | A-posteriori error estimates for the finite element method[END_REF] is proposed and is proved to satisfy Property i) with C = 1 as well as Property ii). In [START_REF] Morin | Local problems on stars: A posteriori error estimators, convergence, and performance[END_REF], an estimator based on the solve of local problems on stars of elements is introduced and it is shown that this estimator satisfies Properties i) and ii) and shows very good accuracy for certain problems such as Poisson problems with an interior layer behavior solution or on a domain with a crack.

A modified version of the Bank-Weiser estimator (based on Neumann problems solved at the level of the cells) that satisfies Property i), Property ii) and Property iv) is introduced in [START_REF] Achchab | Estimateur hiérarchique robuste pour un problème de perturbation singulière[END_REF][START_REF] Achchab | Some remarks about the hierarchicala posteriori error estimate[END_REF] for singularly perturbed reaction-diffusion equations. More recently, novel a posteriori error estimators based on pairs of local Poisson equations with Neumann boundary conditions have been introduced for mixed approximations of the Stokes [START_REF] Liao | A simple yet effective a posteriori estimator for classical mixed approximation of Stokes equations[END_REF] and linear elasticity equations [START_REF] Khan | Robust a posteriori error estimators for mixed approximation of nearly incompressible elasticity[END_REF]. These two estimators are proved to satisfy Property i) and Property ii). In addition, the estimator in [START_REF] Khan | Robust a posteriori error estimators for mixed approximation of nearly incompressible elasticity[END_REF] satisfies Property iv) in the incompressible limit.

In [START_REF] Chen | A PDE approach to fractional diffusion: A posteriori error analysis[END_REF] the implicit error estimator originally introduced in [START_REF] Morin | Local problems on stars: A posteriori error estimators, convergence, and performance[END_REF] is adapted to fractional equations based on the spectral fractional Laplacian operator. This novel estimator is proved to satisfy Property ii).

The Bank-Weiser estimator

In this work, we are particularly interested in the third estimator introduced in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF]. In the rest of this dissertation, we refer to this estimator as the Bank-Weiser estimator. For a cell T , the local contribution of the Bank-Weiser estimator is defined from a hierarchy of Lagrange finite element spaces W T W T . The choice of the polynomial degrees of W T and W T have a strong influence on the properties of the estimator, as well as on what we are able to prove.

xxi

We summarize the state-of-the-art for the Bank-Weiser estimator with a list of properties that have been tackled in the literature:

• Property i) has been established in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF], under a variant of the saturation assumption eq. ( 15),

• a global version of Property ii) has been proved in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF].

• Properties i) and ii) are proved in [START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF] for a particular version of the Bank-Weiser estimator. Hereafter we refer to this estimator as Verfürth's estimator for clarity. As we will see, Verfürth's estimator belongs to a family of Bank-Weiser estimators for which Properties i) and ii) can be established in many cases, without the use of the saturation assumption. The main ingredient being restrictions on the space W ω and in particular the fact that W ω must be of high enough degree.

• In [START_REF] Durán | On the asymptotic exactness of Bank-Weiser's estimator[END_REF], Property iii) is shown for the Bank-Weiser estimator when u h is piecewise linear and on structured triangular meshes, but not in general.

• The Bank-Weiser estimator is shown to satisfy Property i) in [START_REF] Nochetto | Removing the saturation assumption in a posteriori error analysis[END_REF], without the need of a saturation assumption but only when u h is piecewise linear and when Ω is twodimensional.

• In [START_REF] Ainsworth | The performance of Bank-Weiser's error estimator for quadrilateral finite elements[END_REF], some particular cases of the Bank-Weiser estimator are studied and Property iii) is proved for rectangular elements and when u h is a piecewise polynomial function of odd-degree but is proved to be false when u h is of even-degree or when T h is made of non-rectangular elements.

• In [START_REF] Nochetto | Pointwise a Posteriori Error Estimates for Elliptic Problems on Highly Graded Meshes[END_REF], Properties i) and ii) are proved for Verfürth's estimator when • is the pointwise norm.

• In [START_REF] Ainsworth | The influence and selection of subspaces for a posteriori error estimators[END_REF], numerical evidences show that the Bank-Weiser estimator does not satisfy Property iv) with respect to the choice of the spaces W T and W T on rectangular meshes and when u h is a piecewise polynomial of arbitrary degree. In other words, in this context the choice of these spaces has an influence on the performance of the estimator.

• In [START_REF] Verfürth | Robust a posteriori error estimators for a singularly perturbed reactiondiffusion equation[END_REF], a proof is given that Verfürth's estimator satisfy Property iv) with respect to the coefficients of singularly perturbed reaction-diffusion equations.

• In [START_REF] Verfürth | A review of a posteriori error estimation techniques for elasticity problems[END_REF] is shown that Verfürth's estimator, adapted to linear elasticity equations, satisfy Properties i) and ii).

• In [START_REF] Verfürth | Robust A Posteriori Error Estimates for Stationary Convection-Diffusion Equations[END_REF], a proof of Property iv) with respect to the coefficients of stationary convectiondiffusion equations is proposed for Verfürth's estimator.

• In [START_REF] Khan | Robust a posteriori error estimators for mixed approximation of nearly incompressible elasticity[END_REF], the authors show that an estimator based on the solves of local Poisson equations, similar to Bank-Weiser estimator and applied to mixed formulations of the linear elasticity equation satisfy Property iv) in the incompressible limit.

INTRODUCTION

• In Chapter 1, we extend the proof of Property i) for the Bank-Weiser estimator from [START_REF] Nochetto | Pointwise a Posteriori Error Estimates for Elliptic Problems on Highly Graded Meshes[END_REF] to three-dimensional problems.

• In Chapter 2, we provide a novel implementation of the Bank-Weiser estimator in the FEniCS software and give numerical evidences of Property v). In addition, we use the Bank-Weiser estimator to steer adaptive refinement methods.

• In Chapter 3, we propose a method to adapt the Bank-Weiser estimator to fractional equations based on the spectral fractional Laplacian operator for the L 2 error estimation.

Adaptive refinement

In finite element methods the physical domain must be meshed. The mesh plays a crucial role in the accuracy and computational cost of the method. A natural question is: Can we generate a mesh such that the method meets the expected accuracy at a minimal computational cost ? Answering this question is a difficult task and it is the main purpose of adaptive refinement methods. These methods are based on the following feedback loop:

• • • -→ SOLVE -→ ESTIMATE -→ MARK -→ REFINE -→ • • •
This loop is in fact a particular case of fig. 1 in the context of finite element methods. The finite element method is part of our model and the feedback loop of adaptive refinement methods improves the mesh of the finite element method and consequently, improves the model.

Initially a coarse mesh is chosen. The problem is solved (SOLVE) at a very low computational cost on this mesh. An a posteriori error estimator is then used (ESTIMATE) to quantify the discretization error. It is important that this estimator is local. Then depending on the values of the local contributions of the estimator, cells, edges or nodes where the mesh should be refined are selected (MARK). Finally, the zones of the mesh that have been marked during the previous step are refined in order to generate a new mesh. This process is iterated until the global value of the a posteriori error estimator reaches a fixed tolerance, for example. Then, we consider that the mesh is accurate and solve the finite element problem one last time.

Parameters of the method

Adaptive refinement methods rise many questions such as: How should we chose the initial mesh ? What estimator should we chose to steer the method ? What stopping criterion should we use ? How to mark and refine the mesh ? Do adaptive finite element methods converge and if yes, at which rate ?

The origin of adaptive finite element methods goes back [START_REF] Babuška | Error Estimates for Adaptive Finite Element Computations[END_REF] and nowadays most of the above questions have been tackled. A detailed discussion on adaptive refinement methods can be found in [START_REF] Nochetto | Theory of adaptive finite element methods: An introduction[END_REF].

xxiii

The choice of the a posteriori error estimator steering the adaptive method is crucial for the performance of the method. In [START_REF] Carstensen | Estimator Competition for poisson Problems[END_REF] a numerical comparison between various estimators shows that, for Poisson problems, the rate of convergence of the method is usually not influenced by this choice while the global error (and consequently the accuracy of the mesh) can be.

Several marking algorithms are available in the literature such as the maximum strategy, marking the regions with the highest local estimators, proposed in [START_REF] Babuška | Error Estimates for Adaptive Finite Element Computations[END_REF], the bulk chasing strategy, finding the smallest set of marked elements such that the sum of the local estimators is higher than a percentage of the global estimator, introduced in [START_REF] Dörfler | A Convergent Adaptive Algorithm for Poisson's Equation[END_REF] or the equidistribution strategy, which aims at the equidistribution of the local contributions of the error across the mesh, in [START_REF] Morin | A Basic Convergence Result for Conforming Adaptive Finite Elements[END_REF].

One of the challenges of the REFINE step is to preserve the regularity of the mesh, avoiding the creation of cells whose the ratio of the diameter over the diameter of their inscribed ball becomes too large. In addition, the refining process must preserve the conformity of the mesh and avoid the creation of hanging nodes, nodes which are not vertices of their neighboring cells. As for marking strategies, there are multiple refinement algorithms [START_REF] Mitchell | A comparison of adaptive refinement techniques for elliptic problems[END_REF][START_REF] Plaza | Local refinement of simplicial grids based on the skeleton[END_REF][START_REF] Rivara | Mesh Refinement Processes Based on the Generalized Bisection of Simplices[END_REF][START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF], one of the most common being the newest vertex bisection introduced in [START_REF] Mitchell | Optimal Multilevel Iterative Methods for Adaptive Grids[END_REF].

Convergence

The convergence of adaptive refinement algorithms is an important matter and is studied since the late 1990s to nowadays. The first convergence result is proved by Dörfler in [START_REF] Dörfler | A Convergent Adaptive Algorithm for Poisson's Equation[END_REF] under the assumption of a sufficiently fine initial mesh and for two-dimensional problems. An overview of the state-of-the-art on the convergence of adaptive mesh refinement methods is given in [START_REF] Carstensen | Axioms of adaptivity with separate marking for data resolution[END_REF][START_REF] Carstensen | Axioms of adaptivity[END_REF][START_REF] Gantner | Plain convergence of adaptive algorithms without exploiting reliability and efficiency[END_REF].

Fractional Laplacians

In this section we first give a brief sight on fractional calculus, then we discuss different numerical methods used to discretize fractional Laplacian operators and finally we present an overview of the literature on error estimation for finite element discretizations of fractional equations, especially focusing on a posteriori error estimation.

Let -∆ be the usual Laplacian operator. For s ∈ (0, 1), we call fractional Laplacian the fractional power of the Laplacian operator (-∆) s . We are particularly interested in the following fractional equation: let f ∈ L 2 (Ω), find u such that

(-∆) s u = f in Ω, (16a)
u = 0 on ∂Ω. ( 16b
)
For the moment, (-∆) s is just a notation, we need to precise what we mean by this. The fractional Laplacian operator can be defined in multiple ways. In this work we will focus

INTRODUCTION

on the fractional operator called spectral fractional Laplacian and defined from the spectrum of the Laplacian operator on Ω associated to homogeneous zero Dirichlet boundary conditions.

Fractional calculus

The first written appearance of a fractional order derivative comes from a letter Guillaume de l'Hôpital wrote to Gottfried Leibniz in 1695 asking what would be the derivative of order 1/2 of a function. In his answer, Leibniz predicted that this notion would lead to a paradox but many useful consequences would be drawn from this paradox. Leibniz was twice right. Unlike the usual integer order derivative, there is no unique definition of the fractional derivative. Many definitions coexist among others we can cite the Riemann-Liouville derivative [START_REF] Hilfer | Fractional Diffusion Based on Riemann-Liouville Fractional Derivatives[END_REF], the Caputo derivative [START_REF] Caputo | Linear Models of Dissipation whose Q is almost Frequency Independent-II[END_REF] and the Grünwald-Letnikov derivative [START_REF] Scherer | The Grünwald-Letnikov method for fractional differential equations[END_REF]. None of these definitions can take over the others, each one of them having its own advantages and defects. For example, the Grünwald-Letnikov derivative is defined for a very narrow class of functions while the Riemann-Liouville derivative of a constant function is not zero [START_REF] Podlubny | Fractional differential equations : an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF].

Fractional calculus is not only interesting as a mathematical concept but also from an application perspective, one of its main interests being its ability to model non-local phenomena. During the last 50 years, the interest in fractional calculus has grown at an exponential rate [START_REF] Machado | Science metrics on fractional calculus development since 1966[END_REF]. Nowadays, fractional calculus has applications in various fields [START_REF] Sun | A new collection of real world applications of fractional calculus in science and engineering[END_REF] such as e.g. quantum mechanics [START_REF] Laskin | Fractional quantum mechanics and Lévy path integrals[END_REF], economics [START_REF] Song | Solution of the Fractional Black-Scholes Option Pricing Model by Finite Difference Method[END_REF], biology [START_REF] Rihan | Numerical Modeling of Fractional-Order Biological Systems[END_REF], geology [START_REF] Caputo | Linear Models of Dissipation whose Q is almost Frequency Independent-II[END_REF], continuum mechanics [START_REF] Sumelka | Non-local Kirchhoff-Love plates in terms of fractional calculus[END_REF], image processing [START_REF] Antil | Spectral Approximation of Fractional PDEs in Image Processing and Phase Field Modeling[END_REF].

In this work we focus on a particular fractional operator: the fractional Laplacian operator. While the standard non-fractional Lapalcian operator is linked to brownian motion and diffusion, the fractional Laplacian operator is related to Lévy stochastic processes with arbitrary long jumps and anomalous diffusion [START_REF] Bonito | Numerical methods for fractional diffusion[END_REF][START_REF] Pozrikidis | The Fractional Laplacian[END_REF].

Fractional Laplacian definitions and discretizations

When defined on the whole space Ω = R n , many definitions of the fractional Laplacian are equivalent [START_REF] Kwaśnicki | Ten equivalent definitions of the fractional laplace operator[END_REF]. This is no longer the case when the domain Ω is bounded. In this case the way the operator is defined has a significant influence on the solution to the fractional Laplacian equation [START_REF] Duo | A comparative study on nonlocal diffusion operators related to the fractional Laplacian[END_REF]. Moreover, these different operators are infinitesimal generators of different stochastic processes with different physical interpretations [START_REF] Lischke | What is the fractional Laplacian? A comparative review with new results[END_REF]. Two definitions of the fractional Laplacian are particularly studied in the context of the numerical analysis of fractional partial differential equations: the integral fractional Laplacian operator also known as Riesz fractional Laplacian operator [START_REF] Acosta | A Fractional Laplace Equation: Regularity of Solutions and Finite Element Approximations[END_REF][START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF][START_REF] Bonito | Numerical approximation of the integral fractional Laplacian[END_REF][START_REF] Lischke | What is the fractional Laplacian? A comparative review with new results[END_REF] and the spectral fractional Laplacian operator [START_REF] Banjai | Tensor FEM for Spectral Fractional Diffusion[END_REF][START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF][START_REF] Ehstand | Numerical continuation for fractional PDEs: sharp teeth and bloated snakes[END_REF][START_REF] Lischke | What is the fractional Laplacian? A comparative review with new results[END_REF].

There are multiple ways to discretize eqs. (16a) and (16b) when (-∆) s is defined as either the integral or spectral fractional Laplacian operator. Each one of these methods has pros and cons and most of the time the choice of the method will depend on what we want to achieve. xxv Walk-on-spheres method (integral and spectral Laplacians) Originally proposed in [START_REF] Muller | Some Continuous Monte Carlo Methods for the Dirichlet Problem[END_REF], it is a Monte-Carlo method used to simulate paths of s-stable stochastic processes. Since this is a Monte-Carlo method, the main drawback is the slow (logarithmic) convergence rate. However, the main advantages of this method are its very good (linear) complexity scaling as the dimension increases and its "embarrassingly parallel" nature [START_REF] Chati | Random walk method for the two-and three-dimensional Laplace, Poisson and Helmholtz's equations[END_REF].

Grünwald finite differences method (integral Laplacian)

A finite differences based method is proposed in [START_REF] Ciesielski | Numerical Solutions of a Boundary Value Problem for the Anomalous Diffusion Equation with the Riesz Fractional Derivative[END_REF] for one-dimensional problems and modified in [START_REF] Meerschaert | Finite difference methods for two-dimensional fractional dispersion equation[END_REF] for two-dimensional problems. This method shares the same drawbacks as usual finite differences methods and is hardly applicable to domains with complicated boundaries.

Direct (adaptive) finite element method (integral Laplacian)

Using the integral expression of the operator it is possible to derive a weak formulation of eqs. (16a) and (16b) suitable for finite element discretization. This will result in a dense linear system where the matrix entries are based on singular integrals evaluations. A combination of a sparse approximation of the dense matrix based on the panel clustering method and a custom quadrature rule used to compute the singular integrals is proposed in [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF][START_REF] Ainsworth | Towards an Efficient Finite Element Method for the Integral Fractional Laplacian on Polygonal Domains[END_REF], leading to a method with optimal convergence rate and quasi-optimal complexity.

Dirichlet-to-Neumann map method (integral and spectral Laplacians)

This method allows to reformulate the fractional problem in Ω into a non-fractional one with mixed boundary conditions on an extended (semi-infinite cylindrical) domain Ω × R. Such extension results are obtained e.g. in [START_REF] Caffarelli | An Extension Problem Related to the Fractional Laplacian[END_REF][START_REF] Stinga | Extension Problem and Harnack's Inequality for Some Fractional Operators[END_REF]. The Dirichlet-to-Neumann map method is used e.g. in [START_REF] Banjai | Tensor FEM for Spectral Fractional Diffusion[END_REF][START_REF] Bonito | Numerical methods for fractional diffusion[END_REF][START_REF] Chen | A PDE approach to fractional diffusion: A posteriori error analysis[END_REF][START_REF] Nochetto | A PDE Approach to Fractional Diffusion in General Domains: A Priori Error Analysis[END_REF], where a priori and a posteriori error analysis is proposed. The main drawback of such approach is the fact that the discretization requires to work in an extended space with an additional dimension, which makes it difficult to apply on three-dimensional problems.

Spectral element method (spectral Laplacian)

This method consists in a direct approximation of the spectrum of the standard Laplacian using a spectral element method. One of the main advantages of this method is that it is applicable no matter the dimension of the problem. However, it requires the use of quadrilateral meshes which are not always suited for adaptive refinement. Such a method is used e.g. in [START_REF] Song | Computing Fractional Laplacians on Complex-Geometry Domains: Algorithms and Simulations[END_REF].

INTRODUCTION

Eigenproblems solving method (spectral Laplacian)

This method consists in the computation of the full spectrum of the discrete standard Laplacian (discretized using a finite element method). This method takes advantage of the long existing literature on eigenproblems solving numerical methods [START_REF] Ruhe | Rational krylov sequence methods for eigenvalue computation[END_REF]. A scalable algorithm is proposed in [START_REF] Carlson | A scalable framework for solving fractional diffusion equations[END_REF] allowing to solve the eigenvalue problems in a parallel manner. These methods are particularly advantageous when a large number of fractional equations, based on the same spectrum, have to be solved.

(Best) Uniform rational approximation method (spectral Laplacian) This method is based on the approximation of the scalar function λ -→ λ -s by rational functions. They allow to reformulate the fractional problem in a family of non-fractional parametric problems that can be solved independently. The main advantage of such method is the reduced number of parametric problems to solve, however the computation of the rational function coefficients is not always a simple task. Best uniform rational approximations are studied in [START_REF] Harizanov | Optimal solvers for linear systems with fractional powers of sparse SPD matrices[END_REF][START_REF] Harizanov | The Best Uniform Rational Approximation: Applications to Solving Equations Involving Fractional powers of Elliptic Operators[END_REF][START_REF] Harizanov | Numerical solution of fractional diffusion-reaction problems based on BURA[END_REF][START_REF] Hofreither | A unified view of some numerical methods for fractional diffusion[END_REF]. An efficient method to compute coefficients of rational functions, based on barycentric rational interpolation is proposed in [START_REF] Hofreither | An algorithm for best rational approximation based on barycentric rational interpolation[END_REF].

Semi-groups method (integral and spectral Laplacians)

This method is based on the reformulation of the fractional problem into a heat equation which can then be solved using standard techniques such as Euler schemes combined with finite element methods. Like the rational approximation method, it reduces the fractional problem to a family of non-fractional parametric problems. However, in this case the parametric problems cannot be solved independently but must be solved sequentially. Such method is proposed in [START_REF] Cusimano | Discretizations of the Spectral Fractional Laplacian on General Domains with Dirichlet, Neumann, and Robin Boundary Conditions[END_REF][START_REF] Cusimano | Numerical approximations for fractional elliptic equations via the method of semigroups[END_REF].

Reduced order basis method (spectral Laplacian)

As for the rational approximation method, this method proposed in [START_REF] Bonito | Reduced basis approximations of the solutions to spectral fractional diffusion problems[END_REF][START_REF] Danczul | A Reduced Basis Method For Fractional Diffusion Operators I[END_REF][START_REF] Danczul | A Reduced Basis Method For Fractional Diffusion Operators II. may[END_REF][START_REF] Dinh | Model reduction for fractional elliptic problems using Kato's formula[END_REF] is based on the solves of independent non-fractional parametric problems in order to approximate the solution to spectral fractional Laplacian problems. The particularity of this method is to use techniques of model order reduction to optimize the number of parametric problems to solve. This method is particularly attractive when multiple fractional Laplacian equations based on the same spectrum have to be solved.

Dunford-Taylor integral method (integral and spectral Laplacians)

This is the method we are particularly interested in in this work. It is based on the discretizations of integral representations of the fractional Laplacian operators using quadrature rules. Originally developed for spectral fractional elliptic equations in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF], it has since been applied to parabolic equations in [START_REF] Bonito | The approximation of parabolic equations involving fractional powers of elliptic operators[END_REF], space-time fractional parabolic equations in [START_REF] Bonito | Numerical Approximation of Space-Time Fractional Parabolic Equations[END_REF], xxvii fractional powers of regularly accretive operators in [START_REF] Bonito | On sinc quadrature approximations of fractional powers of regularly accretive operators[END_REF] and the integral fractional Laplacian in [START_REF] Bonito | Numerical approximation of the integral fractional Laplacian[END_REF]. This method can be seen as a particular rational approximation method and consists in the reformulation of the fractional problem into a family of non-fractional parametric problems that can be solved independently. We give more details on this method in Chapters 3 and 4.

It is interesting to notice that, for the spectral fractional Laplacian operator, the Dirichletto-Neumann map method, the (best) uniform rational approximation method and the Dunford-Taylor method can all be viewed as rational approximation methods. The study [START_REF] Hofreither | A unified view of some numerical methods for fractional diffusion[END_REF] proposes a unified view on these three methods as particular cases of rational approximation methods as well as numerical comparison of them.

Most of the methods mentioned above are semi-discretization methods and need to be completed to obtain full discretizations. In this work, we focus on a Dunford-Taylor integral method completed by a finite element method. We focus particularly on the discretization error induced by the finite element method in this context.

Error estimation for the fractional Laplacian discretized with finite elements

In this section we give an overview of the literature on error estimation of finite element discretization of fractional equations based on the Laplacian operator, with an emphasis on a posteriori error estimation.

A priori error estimation

The first results on the a priori error estimation of fractional equations finite element discretization date back some fifteen years in [START_REF] Ervin | Variational formulation for the stationary fractional advection dispersion equation[END_REF][START_REF] Ervin | Variational solution of fractional advection dispersion equations on bounded domains in Rd[END_REF][START_REF] Fix | Least squares finite-element solution of a fractional order twopoint boundary value problem[END_REF]. These works focused on fractional advection-diffusion equations. A priori error estimation of the discretization of the integral fractional Laplacian goes back to the a priori estimation results of finite element discretizations of peridynamics equations, see e.g. the review paper [START_REF] D'elia | A Cookbook for Finite Element Methods for Nonlocal Problems, Including Quadrature Rule Choices and the Use of Approximate Balls[END_REF].

A priori error estimation results on the integral fractional Laplacian, seen as a particular case of some non-local operator are proposed in [START_REF] Elia | The fractional Laplacian operator on bounded domains as a special case of the nonlocal diffusion operator[END_REF]. The regularity of the solution to fractional equations based on the integral Laplacian operator is studied in [START_REF] Acosta | A Fractional Laplace Equation: Regularity of Solutions and Finite Element Approximations[END_REF] and optimal convergence rates for the (direct) linear finite element approximation, but only for some values of the fractional power, are derived. The integral Laplacian is discretized using a combination of the Dirichlet-to-Neumann map method (see [START_REF] Caffarelli | An Extension Problem Related to the Fractional Laplacian[END_REF]) and a finite element method in [START_REF] Nochetto | A PDE Approach to Fractional Diffusion in General Domains: A Priori Error Analysis[END_REF] and an a priori error estimation study is proposed. The exponential convergence of the semi-discretization Dirichlet-to-Neumann map method as well as the convergence of the finite element discretization for both quasi-uniform and graded meshes are derived.

At the same time, a priori error estimation results were established in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] for the spectral fractional Laplacian, discretized using the combination of a Dunford-Taylor method and a finite element method. The exponential convergence rate of the Dunford-Taylor method is proved (in the L 2 norm) and the quasi-optimal convergence rate of the finite element method (in the L 2 norm as well) under elliptic regularity assumptions (due to the use of

INTRODUCTION

the Aubin-Nitsche argument). These results were further extended to parabolic equations involving a fractional operator in space [START_REF] Bonito | The approximation of parabolic equations involving fractional powers of elliptic operators[END_REF], space-time fractional parabolic equations [START_REF] Bonito | Numerical Approximation of Space-Time Fractional Parabolic Equations[END_REF] (extension of the semi-discretization Dunford-Taylor method convergence result), steadystate equations but involving a larger class of operators [START_REF] Bonito | On sinc quadrature approximations of fractional powers of regularly accretive operators[END_REF] (and, especially, extending the convergence result to fractional norms) and to the integral Laplacian [START_REF] Bonito | Numerical approximation of the integral fractional Laplacian[END_REF].

The a priori error estimation result from [START_REF] Acosta | A Fractional Laplace Equation: Regularity of Solutions and Finite Element Approximations[END_REF] is improved in [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF][START_REF] Ainsworth | Towards an Efficient Finite Element Method for the Integral Fractional Laplacian on Polygonal Domains[END_REF] and holds for any value of the fractional power in (0, 1), as long as the solution has sufficient interior regularity.

The convergence of spectral fractional Laplacian equations discretizations based on the semigroup method is established in [START_REF] Cusimano | Discretizations of the Spectral Fractional Laplacian on General Domains with Dirichlet, Neumann, and Robin Boundary Conditions[END_REF][START_REF] Cusimano | Numerical approximations for fractional elliptic equations via the method of semigroups[END_REF], for finite elements of arbitrary order and for various homogeneous boundary conditions (Dirichlet, Neumann and Robin).

Using lifting methods, representations of the spectral Laplacian operator associated to non-homogeneous boundary conditions are derived in [START_REF] Antil | Fractional operators with inhomogeneous boundary conditions: analysis, control, and discretization[END_REF] and a priori error estimation results for the Dirichlet-to-Neumann map method combined to a finite element discretization are obtained for problems with non-homogeneous Dirichlet and non-homogeneous Neumann boundary conditions respectively.

The exponential convergence of reduced basis methods applied to the spectral fractional Laplacian is proved in [START_REF] Bonito | Reduced basis approximations of the solutions to spectral fractional diffusion problems[END_REF] and [START_REF] Danczul | A Reduced Basis Method For Fractional Diffusion Operators I[END_REF][START_REF] Danczul | A Reduced Basis Method For Fractional Diffusion Operators II. may[END_REF].

A posteriori error estimation

The first a posteriori error estimator applied to fractional partial differential equations is proposed in [START_REF] Nochetto | A posteriori error analysis for a class of integral equations and variational inequalities[END_REF]. This estimator belongs to the family of explicit residual estimators and is derived for a larger class of integro-differential equations, the integral fractional Lapalcian is part of. The error estimation method is based on bounds of the dual fractional norm by L p norms, allowing to derive an explicit residual estimator based on these norms. Property i) is proved for this estimator and one-dimensional numerical results with efficiency indices are provided.

The work [START_REF] Chen | A PDE approach to fractional diffusion: A posteriori error analysis[END_REF] proposes an a posteriori error estimator for the spectral fractional Laplacian operator. In this work, the fractional operator is discretized using a combination of a Dirichlet-to-Neumann map method (the Caffarelli-Sylvestre extension method [START_REF] Caffarelli | An Extension Problem Related to the Fractional Laplacian[END_REF]) and a finite element method.

Authors of [START_REF] Chen | A PDE approach to fractional diffusion: A posteriori error analysis[END_REF] first claim that explicit residual error estimators cannot be derived in this framework, due to the lack of integration by parts formula (see Section 4.1). Thus, they derive an implicit residual estimator based on local problems on stars of cells (a star associated to a vertex is the set of all cells sharing this vertex), adapting the estimator introduced in [START_REF] Morin | Local problems on stars: A posteriori error estimators, convergence, and performance[END_REF] to the fractional framework. In particular, due to the Dirichlet-to-Neumann map extending the domain with an additional dimension, the local problems are in fact solved on cylindrical stars. One of the main difficulties in this approach is the anisotropy of the mesh in the extended direction. This difficulty is avoided using a geometric condition on the meshes, forcing the cells size in the physical domain direction and size in the extended direction to be of comparable lengths.

Properties i) and ii) are proved for this estimator. However, Property i) stands under the assumption of an unproved conjecture (see Conjecture 5.28 and Remark 5.39). The xxix numerical results provided in the study suggest that this conjecture is in fact true but, to our knowledge, its mathematical justification remains an open question.

It is important to notice that the resulting a posteriori error estimator can only handle the error induced by the finite element discretization and not the error coming from the truncation of the domain of the Caffarelli-Sylvestre extension. In fact, the error due to this truncation decreases exponentially and thus, does not constitutes a bottleneck for the method. However, the only upper bound on this truncation error derived in [START_REF] Chen | A PDE approach to fractional diffusion: A posteriori error analysis[END_REF] depends on a negative fractional Sobolev norm of the data of the problem and is thus not straightforwardly computable. Optimally, an a posteriori error estimator for the truncation error would be required in order to quantify and balance the two sources of error.

Finally, an adaptive refinement method, based on this estimator, the bulk chasing marking strategy and the newest vertex bisection refinement method is designed and numerical results show its quasi-optimal convergence rate.

An explicit residual a posteriori error estimator is proposed in [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF][START_REF] Ainsworth | Towards an Efficient Finite Element Method for the Integral Fractional Laplacian on Polygonal Domains[END_REF]. This estimator is applied to the integral fractional Laplacian and is based on an estimator used in boundary element methods from [START_REF] Faermann | Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary elements methods. Part I. The two-dimensional case[END_REF][START_REF] Faermann | Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary element methods. Part II. The three-dimensional case[END_REF] which is itself derived from the residual estimator introduced in [START_REF] Babuška | A-posteriori error estimates for the finite element method[END_REF]. Property i) and a global version of Property ii) are established in [START_REF] Faermann | Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary elements methods. Part I. The two-dimensional case[END_REF] in the two-dimensional case and in [START_REF] Faermann | Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary element methods. Part II. The three-dimensional case[END_REF] for the three-dimensional case. In addition, Property iv) with respect to the finite element/boundary element polynomial degree is proved.

In [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF][START_REF] Ainsworth | Towards an Efficient Finite Element Method for the Integral Fractional Laplacian on Polygonal Domains[END_REF], the integral fractional Laplacian is discretized using a direct finite element method, leading to a dense stiffness matrix. This stiffness matrix is then approximated by sparse matrices via a cluster method borrowed from the boundary element methods literature. An adaptive refinement method is proposed, based on the residual estimator. One and two-dimensional numerical results are provided showing optimal rates of convergence A gradient recovery error estimation method is proposed in [START_REF] Zhao | Adaptive finite element method for fractional differential equations using hierarchical matrices[END_REF] for one-dimensional fractional differential equations based on the Riesz fractional derivative of order in (1, 2) (which can be then considered as a fractional Laplacian operator). The discretization of this operator by a direct finite element method leads to a dense stiffness matrix which is then approximated by sparse H-matrices [START_REF] Hackbusch | A Sparse Matrix Arithmetic Based on H -Matrices. Part I: Introduction to H -Matrices[END_REF]. The error induced by the approximation of the dense matrix by sparse H-matrices is bounded by a computable quantity that decreases exponentially fast with respect to the approximation parameter.

The a posteriori error estimator is used in an adaptive refinement algorithm based on the bulk chasing marking strategy and a bisection refinement strategy (only one-dimensional problems are considered). Numerical experiments provide evidences of the efficiency of the method on one-dimensional problems as well as the fact that the error estimator satisfy Property v). It is interesting to notice that the errors of the finite element discretization and of the H-matrices method are both controlled and consequently can be balanced.

Fractional diffusion in poro-elasticity

The "real-world" application of this thesis is not concerned with fractional Laplacian equations but with fractional equations based on Caputo's time derivative applied to the anomalous diffusion that occurs in human meniscal tissues.

INTRODUCTION

Caputo's fractional derivative

As we mentioned earlier, the fractional derivative can be defined in many different (and nonequivalent) ways. The Caputo's fractional derivative has been introduced in [START_REF] Caputo | Linear Models of Dissipation whose Q is almost Frequency Independent-II[END_REF] where it is used in dissipative elastodynamics modeling. This derivative has some advantages over other definitions which make it particularly well-suited to fractional differential equations.

The main advantage is the fact that the initial conditions of fractional differential equations based on the Caputo's derivative are conditions on the integer-order derivatives of the solution(s) [START_REF] Podlubny | Fractional differential equations : an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF]. This is particularly useful since it eases the physical interpretation of the initial conditions. Another advantage of Caputo's derivative is that the derivative of a constant function is zero, which is not the case e.g. for the Riemann-Liouville fractional derivative [START_REF] Podlubny | Fractional differential equations : an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF]. We can also mention the fact that, unlike the integer-order derivatives, fractional derivatives are not commutative in general. In other words, if

D α is a fractional derivative of order α ∈ R, then if α = β ∈ R, D α D β f = D β D α f = D α+β f
. However, the conditions of commutativity are less restrictive for the Caputo's derivative compared to e.g. the Riemann-Liouville derivative [START_REF] Podlubny | Fractional differential equations : an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF].

Caputo's derivative found many applications in physics and engineering such as elastodynamics [START_REF] Caputo | Linear Models of Dissipation whose Q is almost Frequency Independent-II[END_REF], fractional relaxation-oscillation and fractional diffusion equations [START_REF] Mainardi | Fractional relaxation-oscillation and fractional diffusion-wave phenomena[END_REF], rheology [START_REF] Mainardi | Creep, relaxation and viscosity properties for basic fractional models in rheology[END_REF] or hydrogeology [START_REF] Iaffaldano | Experimental and theoretical memory diffusion of water in sand[END_REF] and more generally in the study of porous media [START_REF] Barrera | A unified modelling and simulation for coupled anomalous transport in porous media and its finite element implementation[END_REF][START_REF] Caputo | Diffusion in porous layers with memory[END_REF].

Fractional poro-elasticity of the meniscus Among the applications of the Caputo's derivative, we are particularly interested in anomalous diffusion in poro-elastic material and more precisely, in human meniscal tissue. The classical model for diffusion in porous media is given by Darcy's law, stating the proportionality between the flux and the pressure gradient [START_REF] Darcy | Les Fontaines publiques de la ville de Dijon[END_REF][START_REF] Fick | On Liquid Diffusion. The London, Edinburgh and Dublin philosophical magazine and journal of science[END_REF]. However, the classical theory fails to render the behavior of media presenting a spatially and temporally varying permeability. This has for example been highlighted with drained sand layers in [START_REF] Iaffaldano | Experimental and theoretical memory diffusion of water in sand[END_REF], where a model of diffusion integrating memory based on Caputo's fractional derivative [START_REF] Caputo | Models of flux in porous media with memory[END_REF] is validated experimentally. Fractional operators, thanks to their inherent non-locality, are well-suited to model memory effects since memory is equivalent to non-locality in time.

This anomalous behavior is physically interpreted as a consequence of the microstructural rearrangement of the pores during the water flow process, some of the pores being closed, trapping pockets of water and reducing the flow rate while some others remaining open increasing the flow rate [START_REF] Barrera | A unified modelling and simulation for coupled anomalous transport in porous media and its finite element implementation[END_REF][START_REF] Iaffaldano | Experimental and theoretical memory diffusion of water in sand[END_REF].

The (human) meniscus is an example of porous material where this type of behavior is suspected to happen. Meniscal tissues are composed of a porous solid matrix filled with fluid. The material these tissues are composed of is complex, with a non-uniform and anisotropic porosity related to graded material properties [START_REF] Maritz | The functionally grading elastic and viscoelastic properties of the body region of the knee meniscus[END_REF]. The behavior of these tissues is still not well understood [START_REF] Fox | The basic science of human knee menisci: structure, composition, and function[END_REF] which prevents us from conceiving implants successfully mimicking the properties of the meniscus [START_REF] Van Kampen | The knee joint in sports medicine[END_REF].

In Chapter 5, we propose to use the fractional diffusion equation from [START_REF] Caputo | Models of flux in porous media with memory[END_REF] along with (linear) poro-elasticity equations in order to model the behavior of meniscal tissue during xxxi uniaxial confined compression tests. In this particular case, the unidimensional system of equations uncouple and is reduced to time-dependent Poisson equations with analytical solutions [START_REF] Barrera | A unified modelling and simulation for coupled anomalous transport in porous media and its finite element implementation[END_REF], even in the fractional case.

We then validate and calibrate the parameters of our model by fitting experimental results of creep and stress-relaxation tests.

Contributions

In Chapter 1 we provide a novel analysis of Property i) for the Bank-Weiser estimator applied to three-dimensional problems.

The Chapter 2 is concerned with a novel implementation of Bank-Weiser type estimators in the FEniCS finite element software. In addition we propose several numerical results from the application of the estimator to diverse Poisson problems, adaptive refinement algorithms, goal-oriented error estimation as well as a linear elasticity problem. We also provide a time-scale study showing the optimal scaling of our method when implemented using parallel computing.

We adapt the Bank-Weiser method to the estimation of the L 2 error induced by the finite element discretization of the spectral fractional Laplacian in Chapter 3. We propose a method to compute this estimator and to integrate it into an adaptive refinement algorithm for this particular fractional problem. We give evidences of the efficiency of this novel method with several numerical results and especially, provide results from a threedimensional fractional problem. In addition we provide an implementation of this method in the FEniCSx finite element software.

A novel study of the convergence of a rational approximation method is proposed in Chapter 4. From this novel study we deduce a new result, extending an existing result from the literature about the convergence of a particular rational approximation method for the solution to spectral fractional Laplacian equations.

Finally, in Chapter 5 we propose the experimental validation of a fractional poro-elastic model to describe the anomalous diffusion in human meniscal tissues. In addition, calibrate the fractional parameter from this model using a comparison of the model predictions and unidimensional confined compression tests.

INTRODUCTION

Chapter 1

Reliability of the Bank-Weiser estimator in dimension three Introduction

The Bank-Weiser error estimator was introduced in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF]. This seminal work contains a proof that the Bank-Weiser estimator is both efficient and reliable -i.e. it is both a lower and an upper bound of the error-without any restriction on the dimension or on the finite elements order. However, the argument for the upper bound was based on a fragile saturation assumption known to be tricky to assert in practice [START_REF] Dörfler | Small data oscillation implies the saturation assumption[END_REF]. The saturation assumption was successfully removed from the upper bound proof in [START_REF] Nochetto | Removing the saturation assumption in a posteriori error analysis[END_REF] in the case of linear finite elements in dimension two, introducing the additional term referred to as the "data oscillation". An extension of this proof to dimension three does not seem immediate although it is mentioned in the text. In particular, the proof uses the fact that Verfurth's bubble functions [START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF] on edges are quadratic polynomial in dimension two, which is no longer the case in dimension three. In this work, we propose a new proof that is valid both in dimensions two and three. In addition, we provide a short numerical study comparing Bank-Weiser and residual error estimators on a three dimensional test case.

Model problem and finite element discretization

Let Ω ⊂ R d , d = 2 or 3, be a bounded domain with polygonal or polyhedral boundary ∂Ω. For any subdomain ω ⊂ Ω (resp. (d -1)-dimensional set ω), we denote by |ω| the d-dimensional (resp. (d -1)-dimensional) measure of ω. On the domain Ω we consider the usual functions spaces L 2 (Ω) and H 1 0 (Ω). The respective usual norms will be denoted

• ω for L 2 (ω), ω ⊂ Ω and ∇• for H 1 0 (Ω).
For the sake of simplicity, we consider the Poisson equation with homogeneous Dirichlet boundary condition

-∆u = f in Ω, u = 0 on ∂Ω. (1.1)
with given f ∈ L 2 (Ω). The weak form of this problem reads: find u in H 1 0 (Ω) such that for any v in H 1 0 (Ω)

Ω ∇u • ∇v = Ω f v. (1.2)
We discretize this problem using Lagrange piecewise linear continuous finite elements. To do so, we introduce a conformal triangulation T on Ω composed of triangles (resp. tetrahedrons) for d = 2 (resp. d = 3) hereafter called cells. We assume that the triangulation T is regular in the sense of h T /ρ T ≤ γ, ∀T ∈ T , where h T is the diameter of a cell T , ρ T the diameter of its inscribed ball, and γ is positive constant fixed once and for all. For a cell T ∈ T and a non-negative integer p, we denote P p (T ) the set of polynomial functions of degree less than p on T and introduce the spaces of discontinuous and continuous Lagrange finite elements of order p:

V p,dG := {v p ∈ P p (T ), ∀T ∈ T , v p = 0 on ∂Ω} , V p := V p,dG ∩ H 1 0 (Ω).

The finite element approximation to problem eq. (1.2) is: find

u 1 ∈ V 1 such that Ω ∇u 1 • ∇v 1 = Ω f v 1 , ∀v 1 ∈ V 1 . (1.3)
We now introduce some more notations needed in what follows. We call facets the edges of cells in T if d = 2 and the faces of cells in T if d = 3. The notion of triangulation edges will also be important in dimension d = 3. We recall that a facet T ∈ T is a triangle in this case and its boundary consists of 3 sides, called edges. The set of all interior facets of T is denoted by F and the set of all the internal vertices of T is denoted by N . Finally, we use the letter C for various constants that depend only on the triangulation regularity parameter γ and are allowed to change from one occurrence to another.

A posteriori error estimators

Let I : V 2,dG → V 1,dG be the cell by cell Lagrange interpolation operator. The first step in the definition of the Bank-Weiser estimator is to introduce the finite element space V bw := ker(I) = v 2 ∈ V 2,dG , I(v 2 ) = 0 . By definition, the functions of V bw are piecewise quadratic polynomials on the triangulation that vanish at the vertices. Let e bw in V bw be the solution to

T ∈T T ∇e bw • ∇v bw = T ∈T T f v bw + F ∈F F J F {v bw } ∀v bw ∈ V bw , (1.4)
where J F := ∂u 1 ∂n is the jump of the normal derivative of u 1 on F and {•} denotes the average across edges. More precisely, the jump is defined as

J F = (∇u 1 | T F,2 -∇u 1 | T F,1 )
•n where T F,1 , T F,2 are the triangulation cells sharing the facet F and n is the unit normal directed from T F,1 to T F,2 . The Bank-Weiser estimator is then defined as

η bw 2 := T ∈T ∇e bw 2 T . (1.5)
We also recall the explicit residual error estimator

η res 2 := T ∈T h 2 T f 2 T + F ∈F h F J F 2 F . (1.6)
and the data oscillation indicator

osc 2 (f ) := T ∈T h 2 T f -f T 2 T , (1.7)
with h T = diam T and f T = 1 |T | T f . The following theorem establishes the equivalence of the two error estimators eqs. (1.5) and (1.6) modulo a data oscillation term. This proves the equivalence of the Bank-Weiser estimator to the true error ∇u -∇u 1 Ω since such an equivalence is well known to hold for the residual estimator [START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF].

Theorem 1. Let u 1 be the solution of eq. (1.3) for a regular triangulation T . Let η bw and η res be the Bank-Weiser and explicit residual estimators defined respectively in eq. (1.5) and eq. (1.6). Let osc(f ) be the oscillation of f defined in eq. (1.7). Then, there exists two constants c and C only depending on the triangulation regularity such that η bw c η res , (1.8a)

η res C η bw + osc(f ) . (1.8b)
Proof. The arguments to prove eq. (1.8a) can be found in [START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF], applied to a slightly different version of the Bank-Weiser estimator. These arguments consist in putting v bw = e bw in eq. (1.5) and noting (by scaling and equivalence of norms) that v bw T Ch T ∇v bw T and v bw F C √ h T ∇v bw T for any v bw ∈ V bw , T ∈ T and F a facet of T . The proof of eq. (1.8b) essentially proceeds in five steps.

Step 1) Subtracting eq. (1.3) from eq. (1.2) and integrating by parts we get

T ∈T T f v 1 + F ∈F F J F v 1 = Ω ∇(u -u 1 ) • ∇v 1 = 0, ∀v 1 ∈ V 1 .
(1.9)

Taking v 2 in V 2 (continuous) and denoting v 1 = I(v 2 ), v 1 belongs to V 1 and can be used in eq. (1.9). In addition, v 2 -I(v 2 ) belongs to V bw and can be used in eq. (1.4) to get, for all

v 2 ∈ V 2 , Ω ∇e bw • ∇ v 2 -I(v 2 ) = T ∈T T f (v 2 -I(v 2 )) + F ∈F F J F (v 2 -I(v 2 )) = T ∈T T f v 2 + F ∈F F J F v 2 , (1.10)
by linearity of the right hand side and eq. (1.9).

Step 2) For any vertex x ∈ X , one can construct ψ x ∈ V 2 such that ψ x (x) = 1, ψ x = 0 outside of the patch ω x of triangulation cells sharing x and

F ψ x = 0 ∀F ∈ F . (1.11)
In dimension d = 3, we can simply take ψ x = ϕ x the shape function of V 2 associated to x and thus vanishing on the edges midpoints. Indeed, denoting by M F the set of edge midpoints on the facet F (a triangle in this case) we recall that the quadrature rule

F v = 1 3 |F | m∈M F v(m)
is exact on polynomials of degree lower than two. In dimension d = 2, we take ψ x = ϕ x -1 4 m∈Mx ϕ m , where ϕ x is again the shape function of V 2 associated to x, M x is the set of midpoints of the edges sharing x, and ϕ m are the shape functions of V 2 associated to these midpoints. Equation (1.11) is checked for this ψ x by applying Simpson's quadrature rule on the edges.

Using eq. (1.11) and the fact that J F is constant over any facet F , eq. (1.10) with

v 2 = ψ x is reduced to ωx ∇e bw • ∇ ψ x -I(ψ x ) = ωx f ψ x = T ∈ωx f T T ψ x + T ∈ωx T (f -f T )ψ x .
Reordering the terms and using both Cauchy-Schwarz and triangle inequalities give

T ∈ωx f T T ψ x ∇e bw ωx ∇(ψ x -I(ψ x ) ωx + T ∈ωx T (f -f T )ψ x .
(1.12)

Using e.g. quadrature rules for any T in ω x we can compute T ψ

x = -1 6 |T | in dimension d = 2 and T ψ x = -1 20 |T | in dimension d = 3.
Applying the Cauchy-Schwarz inequality twice we also get

T ∈ωx T (f -f T )ψ x T ∈ωx f -f T 2 T 1/2 ψ x ωx .
Moreover, continuity of id -I, Poincaré's inequality as well as a scaling argument give

∇ ψ x -I(ψ x ) ωx C h x |ω x | and ψ x ωx C |ω x |,
where h x is the size of the longest edge in ω x . Then, we finally get

T ∈ωx |T |f T C   1 h x ∇e bw ωx + T ∈ωx f -f T 2 T 1/2   |ω x |. (1.13)
Step 3) Now, for any facet F ∈ F and any cell T ∈ T such that F ⊂ ∂T , one can construct ψ F,T ∈ V bw such that ψ F,T = 1 at all the edge midpoints on F (one midpoint if d = 2 and 3 midpoints if d = 3), ψ F,T = 0 outside of T and

F ψ F,T = 0 ∀F ∈ F , F = F. (1.14)
In dimension d = 2, we take ψ F,T as the usual bubble function associated to the facet F setting ψ F,T = 0 on all the facets of T other than F . In dimension d = 3, we put ψ F,T (m) = -1 2 if m is the midpoint of any edge of T that does not belong to F . Using once again the quadrature rule on triangles F ψ F,T = 1 3 |F | m∈M F ψ F,T (m), we can check that this construction does the job. Now we consider any facet F ∈ F, denote T F,1 , T F,2 the two adjacent triangulation cells and take v bw = ψ F,T F,1 -ψ F,T F,2 in eq. (1.4). The integral of the average {v bw } then vanishes on all the facets and we get

T F,1 ∇e bw • ∇ψ F,T F,1 - T F,2 ∇e bw • ∇ψ F,T F,2 = T F,1 f ψ F,T F,1 - T F,2 f ψ F,T F,2 . (1.15)
Introducing the average of f on cells and reordering the terms give By the precedent quadrature computations, eq. (1.16) and Cauchy-Schwarz inequality we get

f T F,1 T F,1 ψ F,T F,1 -f T F,2 T F,2 ψ F,T F,2 = T F,1 ∇e bw • ∇ψ F,T F,1 - T F,2 ∇e bw • ∇ψ F,T F,2 + T F,1 f T F,1 -f ψ F,T F,1 - T F,2 f T F,2 -f ψ F,T F,2 .
||T F,1 |f T F,1 -|T F,2 |f T F,2 | C    1 h F ∇e bw ω F +   T ∈ω F f -f T 2 T   1/2    |ω F | (1.17) with ω F = T F,1 ∪ T F,2
. Now, if we denote by F x the set of all facets having the vertex x in common and if we consider the finite dimensional vectorial space

E x := {(a T ) T ∈ωx } = R #ωx , the following applications n 1 (a) := T ∈ωx |a T | and n 2 (a) := T ∈ωx a T + F ∈Fx |a T F,1 - a T F,2
|, define norms on E x . Then, using norm equivalence in finite dimension as well as the regularity of the triangulation, we prove the existence of a constant C only depending on triangulation regularity such that

T ∈ωx |T ||f T | C T ∈ωx |T |f T + F ∈Fx |T F,1 |f T F,1 -|T F,2 |f T F,2 .
(1.18)

Step 4) We can now bound the right-hand side of eq. (1.18) using eq. (1.13) and eq. (1.17) to get, for any node x ∈ X ,

T ∈ωx |T ||f T | C   1 h x ∇e bw ωx + T ∈ωx f -f T 2 T 1/2   |ω x | . (1.19)
Taking the square of eq. (1.19) and using triangulation regularity, the fact that f T is constant over the cells and convexity of the square yields

T ∈ωx h 2 T f T 2 T C ∇e bw 2 ωx + T ∈ωx h 2 T f -f T 2 T
.

Summing this over all the vertices and applying once more the triangle inequality as well as triangulation regularity leads to

T ∈T h 2 T f 2 T C η bw 2 + osc 2 (f ) . (1.20)
Step 5) It remains to bound the edge term of the residual estimator eq. (1.6). To this end, we use the functions ψ F,T ∈ V bw again. For any facet F ∈ F , take

v bw = ψ F,T F,1 + ψ F,T F,2 in eq. (1.4) to get J F F v bw T ∈ω F T ∇e bw • ∇v bw + T ∈ω F T f v bw .
Then, by the same quadrature rules as before we have for

T = T F,1 , T F,2 , F ψ F,T = 2 3 |F | in dimension d = 2 and F ψ F,T = |F | in dimension d = 3.
In addition, Cauchy-Schwarz inequality as well as eq. (1.16) give

|F ||J F | C 1 h F ∇e bw ω F + f ω F |ω F | so that h F J F 2 F C( ∇e bw ω F + h 2 F f 2 ω F ). (1.21)
Summing this over all the facets and combining with eq. (1.20) we get eq. (1.8b).

Numerical results

We consider a three-dimensional domain Ω with a L-shaped polyhedral boundary, Ω := (-0.5, 0.5

) 3 \ ([0, -0.5] × [0, -0.5] × [-0.5, 0.5]
). We solve eq. (1.1) with f chosen in order to get the following analytical solution defined on Ω and given, in cylindrical coordinates, by u(r, θ, z) = ϕ(r, θ, z)r 2/3 sin 2θ 3 , where ϕ is a polynomial cut-off function defined (in Cartesian coordinates) by ϕ(x, y, z) = 0.25 -6 (0.25 -x 2 ) 2 (0.25 -y 2 ) 2 (0.25 -z 2 ) 2 . This solution belongs to H 5/3-ε (Ω) for any ε > 0 and its gradient admits a singularity along the re-entrant edge [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF].

We consider here two adaptive refinement algorithms respectively driven by the Bank-Weiser estimator and the residual estimator. Each of these algorithms works as follow:

1. The primal problem eq. ( 1.3) is discretized using piecewise linear Lagrange finite elements.

2. The error is measured using implementations of the estimator (either Bank-Weiser or residual) in the FEniCS Project [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF].

3. The triangulation is marked according to the local contributions of the estimator and using the Dörfler marking strategy. The marking strategy consists in finding the smallest subset M in T such that, T ∈M η 2 T θ 2 η 2 where η = η bw or η res and (η T ) T ∈T are the local contributions of the estimator, for a given parameter θ ∈ (0, 1). Here, we chose the value θ = 0.5 [START_REF] Cascon | Quasi-Optimal Convergence Rate for an Adaptive Finite Element Method[END_REF] in both cases. [START_REF] Aceto | Rational Krylov methods for functions of matrices with applications to fractional partial differential equations[END_REF]. Finally, we refine the triangulation using the Plaza-Carey algorithm present in FEn-iCS [START_REF] Plaza | Local refinement of simplicial grids based on the skeleton[END_REF]. Details on the implementations can be found in [START_REF] Bulle | FEniCS Error Estimation (FEniCS-EE)[END_REF][START_REF] Bulle | An implementation of the Bank-Weiser error estimator in the FEniCS Project finite element software[END_REF][START_REF] Bulle | Hierarchical a posteriori error estimation of Bank-Weiser type in the FEniCS Project[END_REF].

On the left hand side of fig. 1.1 we can see the final triangulation obtained with the Bank-Weiser driven algorithm after four refinement steps. We notice that, as expected, strong refinement occurs near the re-entrant corner edge around the origin. The choice of the a posteriori error estimator does not have a strong influence on the mesh hierarchy.

On the right of fig. 1.1, the convergence curves of the estimator and the corresponding exact error are plotted for both refinement algorithms. The approximated rates of convergence (estimated with least squares) are, in the case of the Bank-Weiser driven adaptive algorithm, -0.33 for the exact error and -0.31 for the Bank-Weiser estimator and, in the case of the residual driven adaptive algorithm, -0.34 for the exact error and -0.33 for the residual estimator.

In the case of the Bank-Weiser driven adaptive algorithm, we have also computed the residual estimator on the same mesh hierarchy in order to compare estimators efficiencies η/ ∇(u -u h ) for η = η bw or η res . Estimators efficiencies on the last refinement step are respectively 1.9 for the Bank-Weiser estimator and 5.9 for the residual estimator. We can notice that the Bank-Weiser estimator is much sharper than the residual estimator but it is not as sharp as for two-dimensional problems (see e.g. [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF]). The fact that the Bank-Weiser estimator is not asymptotically exact for non-structured meshes is known and was proved in [START_REF] Durán | On the asymptotic exactness of Bank-Weiser's estimator[END_REF].

Introduction

A posteriori error estimation [START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF] is the de facto tool for assessing the discretization error of finite element method (FEM) simulations, and iteratively reducing that error using adaptive mesh refinement strategies [START_REF] Nochetto | Theory of adaptive finite element methods: An introduction[END_REF].

This paper is concerned with the description and justification of an implementation of an error estimator introduced in the seminal paper of Bank and Weiser [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF]Section 6]. In that paper an error estimate was derived involving the solution of local Neumann problems on a special finite element built on nested or hierarchical spaces. Despite its excellent performance and low computational cost, this estimator has seen relatively sparse use in practical computational problems. The overarching goal of this contribution is to provide access to an efficient, generic and extensible implementation of Bank-Weiser type estimators in a modern and widely used finite element software, specifically, the FEniCS Project [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF].

Background

The literature on a posteriori error estimation and adaptive finite element methods is vast, so we focus on articles on practical software implementations of adaptive finite element methods and comparative performance studies.

The T-IFISS [START_REF] Bespalov | T-IFISS: a toolbox for adaptive FEM computation[END_REF] software package, based on the existing IFISS [START_REF] Elman | IFISS: A Computational Laboratory for Investigating Incompressible Flow Problems[END_REF] package, is a finite element software written in MATLAB/Octave with a focus on a posteriori error estimation and adaptive finite element methods. Recently [START_REF] Bespalov | Goal-oriented error estimation and adaptivity for elliptic PDEs with parametric or uncertain inputs[END_REF], T-IFSS has been extended to solve adaptive stochastic Galerkin finite element methods. The stated emphasis of T-IFISS [START_REF] Bespalov | T-IFISS: a toolbox for adaptive FEM computation[END_REF] is on being a laboratory for experimentation and exploration, and also to enable the rapid prototyping and testing of new adaptive finite element methods. A number of estimation and marking strategies are implemented in T-IFISS, although not the Bank-Weiser estimator we consider in this paper. T-IFISS only works for two-dimensional problems and it was never intended to be a high-performance code suitable for large-scale computations e.g. high-performance computing systems using the Message Passing Interface (MPI).

The PLTMG package [START_REF] Bank | PLTMG: A Software Package for Solving Elliptic Partial Differential Equations: Users' Guide 8[END_REF] is one of the oldest open finite element softwares for solving elliptic problems that is still under active maintenance, and includes many advanced features such as hp-adaptive refinement, a posteriori error estimation, domain decomposition and multigrid preconditioning. The a posteriori error estimation is based on a superconvergent patch recovery estimation technique introduced in [START_REF] Bank | Superconvergent Derivative Recovery for Lagrange Triangular Elements of Degree p on Unstructured Grids[END_REF]. PLTMG only works in two dimensions and is naturally limited from a usability perspective due to the programming tools available at its inception (Fortran and ANSI C).

In [START_REF] Funken | Efficient implementation of adaptive P1-FEM in Matlab[END_REF] an adaptive first-order polynomial finite element method was implemented in a code called p1afem using MATLAB. The primary goal was to show how the basic finite element algorithm could be implemented efficiently using MATLAB's vectorization capabilities. A standard residual estimator [START_REF] Babuška | Feedback and adaptive finite element solution of onedimensional boundary value problems[END_REF] is used to drive an adaptive mesh refinement algorithm. Again, like T-IFISS, p1afem only works in two dimensions.

In [START_REF] Rognes | Automated Goal-Oriented Error Control I: Stationary Variational Problems[END_REF] a novel methodology for automatically deriving adaptive finite element methods from the high-level specification of the goal functional and (potentially non-linear) residual equation was implemented in the FEniCS Project. The emphasis of the paper [START_REF] Rognes | Automated Goal-Oriented Error Control I: Stationary Variational Problems[END_REF], in contrast with the T-IFISS toolbox [START_REF] Bespalov | T-IFISS: a toolbox for adaptive FEM computation[END_REF], is on the automatic construction of goal-oriented adaptive finite element methods, without much knowledge required on the part of the user. The implicit residual problems are automatically localized using bubble functions living on the interior and facets of the cell, and the dual problem [START_REF] Giles | Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality[END_REF] is derived and solved automatically on the same finite element space as the primal problem, before being extrapolated to a higher-order finite element space using a patch-wise extrapolation operator. In practice the automatically derived estimators seem to be able to effectively drive adaptive mesh refinement for a range of different PDEs.

Explicit residual estimators are also commonly employed by users of high-level finite element software packages as they can usually be expressed straightforwardly in a high-level form language, e.g. [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF][START_REF] Prud'homme | Feel++ : A computational framework for Galerkin Methods and Advanced Numerical Methods[END_REF]. For example, [START_REF] Hale | Simple and extensible plate and shell finite element models through automatic code generation tools[END_REF] used the FEniCS Project to implement an explicit residual error estimator for the Reissner-Mindlin plate problem from [START_REF] Beirão Da Veiga | A-priori and aposteriori error analysis for a family of Reissner-Mindlin plate elements[END_REF]. The authors of [START_REF] Duprez | Quantifying discretization errors for soft tissue simulation in computer assisted surgery: A preliminary study[END_REF] used the FEniCS Project to implement an explicit residual estimator for elasticity problems within a dual-weighted residual framework. The dual problem is solved on a higher-order finite element space in order to ensure that the weighting by the dual residual solution does not vanish [START_REF] Rognes | Automated Goal-Oriented Error Control I: Stationary Variational Problems[END_REF]. In [START_REF] Houston | Automatic Symbolic Computation for Discontinuous Galerkin Finite Element Methods[END_REF] the authors use an explicit dual-weighted residual strategy for adaptive mesh refinement of discontinuous Galerkin finite element methods. In addition, as the name suggests, they can be explicitly computed as they involve only functions of the known finite element solution and the problem data.

In the present work, aside of the Bank-Weiser estimator we will consider an explicit residual estimator [START_REF] Babuška | A-posteriori error estimates for the finite element method[END_REF] named residual estimator in the following, a flux reconstruction based on averaging technique estimator [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF], referred to as Zienkiewicz-Zhu estimator, and a variant of the Bank-Weiser estimator introduced in [START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF] and referred to as the bubble Bank-Weiser estimator. The residual estimator was proved to be both reliable and (locally) efficient in [START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF] for any finite element order and in any dimension. The proof of reliability and (local) efficiency of Zienkiewicz-Zhu estimator has been derived in [START_REF] Rodríguez | Some remarks on Zienkiewicz-Zhu estimator[END_REF], for linear finite elements in dimension two and generalized to any averaging technique in any dimension in [START_REF] Carstensen | Each averaging technique yields reliable a posteriori error control in FEM on unstructured grids. Part I: Low order conforming, nonconforming, and mixed FEM[END_REF] and any finite element order in [START_REF] Bartels | Each averaging technique yields reliable a posteriori error control in FEM on unstructured grids. Part II: Higher order FEM[END_REF]. The bubble Bank-Weiser estimator was proved to be reliable and locally efficient in [START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF] for any dimension and any finite element order.

A proof of the equivalence between the Bank-Weiser estimator and the exact error was derived in the original paper [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF]. However, this proof requires a saturation assumption [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF][START_REF] Dörfler | Small data oscillation implies the saturation assumption[END_REF][START_REF] Nochetto | Removing the saturation assumption in a posteriori error analysis[END_REF] asking for the best approximation with higher order finite elements to be strictly smaller than that of lower order elements and which is known to be tricky to assert in practice. Some progress has been made in [START_REF] Nochetto | Removing the saturation assumption in a posteriori error analysis[END_REF] removing the saturation assumption from the analysis. However, this progress was made at the price of restricting the framework to linear polynomial finite elements and dimension two only. The equivalence proof between Bank-Weiser and residual estimators have been improved by the authors in [START_REF] Bulle | Removing the saturation assumption in Bank-Weiser error estimator analysis in dimension three[END_REF] where it was extended to dimension three.

Contribution

We show how robust and cheap hierarchical error estimation strategies can be implemented in a high-level finite element framework, e.g. the FEniCS Project [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF], Firedrake [START_REF] Gibson | Slate: extending Firedrake's domain-specific abstraction to hybridized solvers for geoscience and beyond[END_REF][START_REF] Rathgeber | Firedrake: Automating the Finite Element Method by Composing Abstractions[END_REF], freefem++ [START_REF] Hecht | New development in freefem++[END_REF], Feel++ [START_REF] Prud'homme | Feel++ : A computational framework for Galerkin Methods and Advanced Numerical Methods[END_REF], GetFEM [START_REF] Renard | GetFEM: Automated FE modeling of multiphysics problems based on a generic weak form language[END_REF] or Concha [START_REF]Complex Flow Simulation Codes based on Highorder and Adaptive methods[END_REF]. Specifically, the contribution of our paper to the existing literature is:

• A generic and efficient implementation of the Bank-Weiser estimator in the open source FEniCS Project finite element software that works for Lagrange finite elements of arbitrary polynomial order and in two and three spatial dimensions. We provide implementations for the popular but legacy DOLFIN finite element solver [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF], and the new DOLFINx solver [START_REF] Habera | FEniCSX: A sustainable future for the FEniCS Project[END_REF]. The two versions are functionally identical, although in terms of overall speed and parallel scaling the DOLFINx version is superior due to underlying architectural improvements. Hence we only show parallel scaling results with this new version. The code is released under an open source (LGPLv3) license [START_REF] Bulle | An implementation of the Bank-Weiser error estimator in the FEniCS Project finite element software[END_REF]. Because the code utilizes the existing automatic code generation capabilities of FEniCS along with a custom finite element assembly routine, the packages are very compact (a few hundred lines of code, plus documentation and demos). Additionally, the estimators are implemented in near mathematical notation using the Unified Form Language, see Sections 2.7.1 and 2.8.1 for code snippets.

• A numerical comparison of the Bank-Weiser estimator with various estimators mentioned earlier. We examine the relative efficiency, and their performance within an adaptive mesh refinement loop on various test problems. Unlike [START_REF] Carstensen | Estimator Competition for poisson Problems[END_REF], we do not aim at running a competition of error estimators but at stressing the potential of the Bank-Weiser estimator since, as the authors of [START_REF] Carstensen | Estimator Competition for poisson Problems[END_REF] point out, a single error estimation strategy is not sufficient to cover the particulars of all possible problems.

• Relying on results in [START_REF] Becker | Weighted Marking for Goal-oriented Adaptive Finite Element Methods[END_REF], we show a goal-oriented adaptive mesh refinement algorithm can be driven by weighted sum of estimators, computed separately on primal and dual problems discretized on the same finite element space. This avoids the extrapolation operation of [START_REF] Rognes | Automated Goal-Oriented Error Control I: Stationary Variational Problems[END_REF] or the need to compute the dual solution in a higherorder finite element space [START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF].

• Using the same basic methodology as for the Poisson problem, we extend our approach to estimating errors in mixed approximation of nearly incompressible elasticity problems. This idea was originally introduced in [START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF] and is still an active research topic, see e.g. [START_REF] Khan | Robust a posteriori error estimators for mixed approximation of nearly incompressible elasticity[END_REF] for a parameter-robust implicit residual estimator for nearlyincompressible elasticity.

Outline

An outline of this paper is as follows; in Section 2.1.4 we outline the main notation and definitions used in this paper. In Sections 2.2 and 2.3 we show the derivation of the primal problem and the Bank-Weiser error estimator. In Section 2.4 we derive a new method for computing the Bank-Weiser estimator and discuss its implementation in FEniCS. In Section 2.5 we discuss the use of the approach for various applications such as goal-oriented adaptive mesh refinement and for mixed approximations of PDEs. Then, in Section 2.6 we show some results on two and three dimensional Poisson test problems as well as on linear elasticity problems, before concluding in Section 2.10.

Notation

In this section we outline the main notations used in the rest of the paper. Let Ω be an bounded open domain of R d (d = 1, 2 or 3), with polygonal/polyhedral boundary denoted by Γ := ∂Ω. We consider Γ = Γ D ∪ Γ N a partition of the boundary. We assume Γ D is of positive measure. We denote by n : Γ → R d the outward unit normal vector along Γ.

Let ω be a subset of Ω. For l ∈ R we denote by H l (ω) the Sobolev space of order l. The space H 0 (ω) = L 2 (ω) is the Lebesgue space of square integrable functions over ω. The space H l (ω) is endowed with the usual inner product (•, •) l,ω and norm • l,ω . We omit the subscript l when l = 0 and subscript ω when ω = Ω. We denote H 1 D (Ω) the subspace of H 1 (Ω) of functions with zero trace on Γ D . We make use of the notation ∂ n v := ∇v • n for the normal derivative of a smooth enough function v. For l ∈ R and for a d-dimensional subset ω of Ω, we also define the following vector fields spaces L 2 (ω) := L 2 (ω)

d and

H l (ω) := H l (ω)
d , with respective inner products defined as their scalar counterparts, replacing the scalar product by the Euclidean inner product or the Frobenius double dot product. The space H 1 D (Ω) is the subspace of H 1 (Ω) of functions with zero trace on Γ D . From now on, the bold font notation will be reserved to vector fields. With these notations at hand we can proceed with the rest of the paper.

Primal problem statement and finite element discretization

-∆u = f in Ω, u = u D on Γ D , ∂ n u = g on Γ N . (2.1)
Problem eq. (2.1) can be written in an equivalent weak form:

Find u ∈ H 1 (Ω) of trace u D on Γ D such that (∇u, ∇v) = (f, v) + (g, v) Γ N , ∀v ∈ H 1 D (Ω). (2.2)
The weak problem eq. ( 2.2) can be discretized using the Lagrange finite element method. We take a mesh T of the domain Ω, consisting of cells T = {T }, facets E = {E} (we call facets the edges in dimension two and the faces in dimension three), and vertices N = {χ}. The mesh T is supposed to be regular in Ciarlet's sense: h T /ρ T γ, ∀T ∈ T , where h T is the diameter of a cell T , ρ T the diameter of its inscribed ball, and γ is a positive constant fixed once and for all. The subset of facets in the interior of the mesh (i.e. those that are not coincident with the boundary Γ) is denoted E I . The subset of facets lying on

Γ D is denoted E D .
The subset of facets lying on Γ N is denoted E N . The subset of facets lying on the boundary of the domain Γ is denoted

E B = E D ∪ E N .
Until the end of this work we assume that the mesh resolves the boundary conditions, in other words for any edge

E ∈ Γ then E ∈ Γ D or E ∈ Γ N . Let n + ∈ R d and n -∈ R d
be the outward unit normals to a given edge as seen by two cells T + and T -incident to a common edge E. If we denote P k (T ) the space of polynomials of order k on a cell T , the continuous Lagrange finite element space of order k on the mesh T is defined by

V k := v k ∈ H 1 (Ω), v k|T ∈ P k (T ) ∀T ∈ T . (2.3)
We denote V k D the finite element space composed of functions of V k vanishing on the boundary Γ D . We consider the finite element problem: Find

u k ∈ V k such that u k = u D,k on Γ D and: (∇u k , ∇v k ) = (f, v k ) + (g, v k ) Γ N , ∀v k ∈ V k D , (2.4)
and where u D,k is a discretization of u D on V k (for example the Laplace interpolation or a L 2 orthogonal projection).

The Bank-Weiser estimator

In this section we derive the general definition of the Bank-Weiser estimator from the equation of the error as it was given in the original paper [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF]. We also give a concrete example of the Bank-Weiser estimator for linear finite elements.

The global error equation

We are interested in estimating the error we commit by approximating the solution u by

u k ∈ V k D .
We define this error by the function e := u -u k and we want to estimate its norm e 1 . The first step towards this will be to derive a new variational problem for which the exact error e is the solution. For a cell T of the mesh, we introduce the interior residual as

r T := (f + ∆u k ) |T , (2.5)
and for an edge E, the edge residual

J E =    0 if E ∈ E D , ∂ n u k E if E ∈ E I , (g -∂ n u k ) |E if E ∈ E N . (2.6)
where the notation v E := v + -v -denotes the jump in the value of the function across an interior facet E ∈ E I . Here, v + and v -denote the values of v on the facet E as seen by the two incident cells T + and T -, respectively. The error function e satisfies what we call the global error equation

(∇e, ∇v) = T ∈T (r T , v) T + E∈E I (J E , v) E + E∈E N (J E , v) E , ∀v ∈ H 1 D (Ω), (2.7)
and e = u D -u k on the Dirichlet boundary Γ D .

The local Bank-Weiser space and the Bank-Weiser estimator

We introduce now local finite element spaces in order to derive the finite element approximation of the error. For a cell T of the mesh we define

V k T,D := v k,T ∈ P k (T ), v k,T = 0 in (Ω \ T ) ∪ (T ∩ Γ D ) , (2.8)
as well as

V k T := {v k,T ∈ P k (T )} .
(2.9)

A key idea in the Bank-Weiser estimator derivation is to introduce an appropriate finite element space for the discretization of error. This non-standard space has two roles. Firstly, for the local problems involving the cells with facets only in the interior of the domain or on the Neumann boundary, it should remove the constant functions, giving a unique solution. Secondly, and as we will notice in Section 2.6, solving the local error equation on the finite element space V k T,D /R does not necessary lead to an accurate estimation of the error. However, in some cases, the estimation of the error can be surprisingly accurate when the space is judiciously chosen. We refer the reader to [START_REF] Ainsworth | The influence and selection of subspaces for a posteriori error estimators[END_REF] for a full discussion.

Before introducing this non-standard space, we need some more notations. Let k + and k -be two non-negative integers such that k + > k -0. Let T be the reference cell fixed once for all (independent from the mesh T ). We denote

L T : V k + T -→ V k + T , Im(L T ) = V k - T , (2.10)
the Lagrange interpolation operator between the local spaces

V k + T and V k - T ⊂ V k + T
. Moreover, for any cell T of the mesh, there exists an affine bijection

S : T -→ T x -→ S( x) =: x (2.11)
mapping T onto T . From the mapping S we deduce another mapping given by

S : V k + T -→ V k + T v(x) -→ S(v)( x) := v(S( x)).
(2.12)

If we denote d + the dimension of V k + T and d -the dimension of V k - T , given B + T := { ϕ 1 , • • • , ϕ d + } the basis of shape functions of V k + T and B + T := {ϕ 1 , • • • , ϕ d + } the basis of V k +
T , we can always find a mapping S (and a mapping S) such that

S(ϕ T,i ) = ϕ T,i , ∀i ∈ {1, • • • , d + }, (2.13)
We choose S and S so. For a given cell T of the mesh, we define the Lagrange interpolation operator on T as follows

L T := S -1 • L T • S. ( 2
G = Id -1 G Id = G.
(2.15)

For a cell T of the mesh, the local Bank-Weiser space V bw T is defined as the null space of L T , in other words

V bw T := ker(L T ) = v bw T ∈ V k + T , L T v bw T = 0 . (2.16)
Similarly, we define

V bw T,D := v bw T ∈ V bw T , v bw T = 0 on T ∩ Γ D .
(2.17)

With these new spaces in hands, we can derive a local discrete counterpart of eq. (2.7) on any cell T : Find e T ∈ V bw T such that:

∇e bw T , ∇v bw T = r T , v bw T + 1 2 E∈∂T J E , v bw T E , ∀v bw T ∈ V bw T,D , (2.18
)

and e T = π bw T (u D -u k ) on Γ D , where π bw T : L 2 (T ) → V bw
T is a proper projection operator (the way this projection is implemented is detailed in Section 2.4.1).

Note, the definition of the edge residual J takes into account the error on the Neumann boundary data approximation. The Dirichlet boundary data approximation has to be incorporated to the linear system during the solve of eq. (2.18), as well will see later. For a detailed discussion on a priori and a posteriori error estimation with inhomogeneous Dirichlet boundary conditions see [START_REF] Aurada | Each H 1/2 -stable projection yields convergence and quasi-optimality of adaptive FEM with inhomogeneous Dirichlet data in R d[END_REF][START_REF] Bartels | Inhomogeneous Dirichlet conditions in a priori and a posteriori finite element error analysis[END_REF].

Finally, on the cell T the local Bank-Weiser estimator η bw,T is defined by

η bw,T := ∇e bw T T , (2.19)
where e T is defined in eq. (2.18) and the global Bank-Weiser estimator by the sum of local estimates

η 2 bw := T ∈T η 2 bw,T . (2.20)
Note, although it is not shown in this study, it is straightforward to generalize the Bank-Weiser estimator for other kind of elliptic operators by changing the energy norm in eq. (2.19) accordingly.

A particular example

If we assume k = 1 (i.e. we solve eq. (2.4) using linear finite elements) one can define the space V bw T from the choice of k + = 2, k -= k = 1. This example was the case considered in the numerical tests of the original paper [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF]. The space V bw T consists of quadratic polynomial functions (in V 2

T ) vanishing at the degrees of freedom of the standard linear finite element functions (in V 1 T ) i.e. the degrees of freedom associated with the vertices of T .

Algorithms and implementation details

The linear system corresponding to eq. (2.18) is not accessible in FEniCS. This prevent us from directly solving the Bank-Weiser equation. We propose to bypass the problem by constructing the linear system corresponding to eq. (2.18) from another linear system derived from finite element spaces that are accessible directly in FEniCS.

Method outline

1. We consider the following singular value decomposition (SVD) of

G G = U ΣV T , (2.21)
where Σ is a diagonal matrix composed of the singular values of G. The columns of the matrix V are singular vectors of G, associated with singular values. The columns associated with singular values zero span the null space of G. We take the submatrix N made of the columns of V spanning the null space of G. Note that, since G does not depend on any cell T , the same property holds for N .

2. We build the matrix A + T and vector b + T of the local linear system corresponding to the following variational formulation in the space V k + T , available in FEniCS:

∇e + T , ∇v + T = r T , v + T + 1 2 E∈∂T J E , v + T E , ∀v + T ∈ V k + T . (2.22)
We integrate the Dirichlet boundary condition directly into A + T and b + T , by considering the vector associated to π + T (u D -u k ), where π + T is the L 2 projection onto V k + T . More precisely, the rows and columns of A + T corresponding to degrees of freedom on the Dirichlet boundary are zeroed and the corresponding diagonal entries are replaced by ones. The entries of b + T corresponding to these degrees of freedom are replaced by the corresponding entries in the vector of π + T (u D -u k ).

3. We construct the matrix A bw T and vector b bw T as follow

A bw T = N T A + T N and b bw T = N T b + T , (2.23)
where A bw T and b bw T are the matrix and vector which allow to recover the bilinear and linear forms of eq. (2.18) in a basis of V bw T .

4. We solve the linear system

A bw T x bw T = b bw T , (2.24)
5. We bring the solution back to V k + T , considering N x bw T , in order to post-process it and compute the local contribution of the Bank-Weiser estimator eq. (2.19).

Compute G = G 1 G 2 SVD of G to obtain V Extract N from V Computation of N Compute A + T and b + T Compute A bw

Computational details

We now give more details specific to our implementation in FEniCS of each one of the above steps.

1. Computation of N . This is the key point of our implementation. The operator L T can be written as follows:

L T : V k + T -→ V k - T -→ V k + T v + -→ G 1 (v + ) -→ G 2 G 1 (v + ) . (2.25)
Then, the matrix G is obtained via the following product

G = G 2 G 1 , (2.26)
where G 1 and G 2 are respectively the matrix in the couple of basis

(B + T , B - T ) of the Lagrange interpolation operator from V k + T to V k - T , denoted G 1 and the matrix in the same couple of basis of the canonical injection of V k - T into V k + T , denoted G 2 .
The matrices G 1 and G 2 can be calculated either using the Finite Element Automatic Tabulator (FIAT) [START_REF] Kirby | Algorithm 839: FIAT, a new paradigm for computing finite element basis functions[END_REF] or, as we choose to do, using the interpolator construction functions of the DOLFIN/x finite element library [START_REF] Logg | DOLFIN: Automated finite element computing[END_REF]. The next step consists in computing the unitary matrix V of right singular vectors of G. This computation is done using the singular value decomposition (SVD) algorithm available in the SciPy library [START_REF] Virtanen | and Contributors. SciPy 1.0-Fundamental Algorithms for Scientific Computing in Python[END_REF]. We can write the matrix V as follows,

V = ξ 0 1 | • • • | ξ 0 d bw | ξ 1 | • • • |ξ d -, (2.27)
where

B bw T := {ξ 0 1 , • • • , ξ 0 d bw } is the set of singular vectors of G corresponding to a zero singular value, spanning V bw T and {ξ 1 , • • • , ξ d -} is spanning the supplementary space.
The matrix N is then chosen as the submatrix of V , keeping only the columns from B bw T :

N := ξ 0 1 | • • • | ξ 0 d bw .
(2.28)

The linear algebra operations needed to form the submatrix N from V are performed using the NumPy library [START_REF] Van Der Walt | The NumPy Array: A Structure for Efficient Numerical Computation[END_REF].

Computation of A +

T and b + T . The eq. (2.22) is expressed directly in the Unified Form Language (UFL) [START_REF] Alnaes | Unified form language: A domain-specific language for weak formulations of partial differential equations[END_REF] and efficient C++ code for calculating the cell local tensors A + T and b + T for a given cell T is then generated using the FEniCS Form Compiler (FFC) [START_REF] Kirby | A compiler for variational forms[END_REF][START_REF] Ølgaard | Automated Code Generation for Discontinuous Galerkin Methods[END_REF]. If the cell T has an edge on a Dirichlet boundary E D , the matrix A + T and vector b + T must be modified in order to enforce the boundary condition.

Computation of A bw

T and b bw T . The matrix A bw T and vector b bw T are constructed using eq. (2.23).

4. Solution of the linear system (2.24). The linear system eq. (2.24) is solved using a partialpivot LU decomposition algorithm from the Eigen dense linear algebra library [149] in DOLFIN and xtensor-blas, which calls LAPACK's dgesv in DOLFINx.

5.

Computation of the Bank-Weiser estimator. Finally, the solution x bw T is sent back to V k + T using N and the norm of the corresponding function, giving the local estimator eq. (2.19) is computed using standard high-level functions already available within FEn-iCS. The global estimator eq. (2.20) is computed using the information of all the local contributions.

Additional remarks

• The custom assembler composed of steps 2.-5. is performed by looping over every cell of the mesh and, by virtue of using the abstractions provided by DOLFINx, works in parallel on distributed memory computers using the Message Passing Interface (MPI) standard. For performance reasons these steps have been written in C++ and wrapped in Python using the pybind11 library so that they are available from the Python interface to DOLFIN/x. In contrast, the first step must only be performed once since the matrix N is the same for every cell of the mesh.

• A posteriori error estimation methods such as the one we are considering here assume that the linear system associated with the primal problem eq. (2.2) is solved exactly. However for performance reasons, here we use PETSc conjugate gradient iterative method. Using inexact solutions can have an influence on the total error but also on the a posteriori error estimator itself. It is a known issue [START_REF] Arioli | Interplay between discretization and algebraic computation in adaptive numerical solutionof elliptic PDE problems[END_REF] and several authors have proposed ways to estimate the algebraic error, see e.g. [START_REF] Anciaux-Sedrakian | Adaptive solution of linear systems of equations based on a posteriori error estimators[END_REF][START_REF] Papež | Estimating and localizing the algebraic and total numerical errors using flux reconstructions[END_REF]. Since algebraic error estimation is beyond the scope of this work, in all our numerical results we set PETSc residual tolerance small enough to neglect this part of the error.

• Because we use the automatic code generation capabilities of FEniCS, our approach can be readily applied to other definitions for the spaces

V k + T and V k -
T , and to vectorial problems like linear elasticity, as we will see in the next section.

• For large problems the storage of the global higher order space V k + can be an issue since it requires a lot of memory space. However we avoid this problem by considering the local higher order spaces V k + T

(and local lower order spaces V k - T ) only.

• In the numerical results section we compare several versions of Bank-Weiser estimator and especially the one we call bubble Bank-Weiser estimator and denote η b T which can be obtained with our method by taking V + T as the space V 2 T + Span{ψ T } (the local space of quadratic functions enriched with the space spanned by the interior bubble function) and

V k - T as V 1
T . The resulting space V bw T is spanned by the interior bubble function and the edges bubbles functions of the cell T .

Applications

In this section we show a number of applications, including adaptive mesh refinement, goal-oriented estimation and extensions to more complex mixed finite element formulations for the nearly-incompressible elasticity problems.

Adaptive mesh refinement

As well as simply providing an estimate of the global and local error, the estimator can be used to drive an adaptive mesh refinement strategies. In the following we compare different refinement strategy all based on the following loop:

• • • -→ SOLVE -→ ESTIMATE -→ MARK -→ REFINE -→ • • •
The loop can be terminated once a given criterion e.g. maximum number of iterations, or global error less than a given tolerance, has been reached. A detailed discussion on adaptive refinement methods can be found in [START_REF] Nochetto | Theory of adaptive finite element methods: An introduction[END_REF]. In the following we expand on the specific algorithms used in our case.

Solve

The weak form eq. (2.2) is discretized using a standard finite element method implemented within FEniCS. The resulting linear systems are solved using the appropriate algorithms available within PETSc [START_REF] Balay | PETSc Users Manual[END_REF], e.g. conjugate gradient method preconditioned with Hypre BoomerAMG [START_REF] Falgout | hypre: A Library of High Performance Preconditioners[END_REF], or direct methods, e.g. MUMPS [START_REF] Amestoy | A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF].

Estimate

The Bank-Weiser estimator η bw is formulated and implemented as described in Section 2.4. The local contributions of the estimator provide an estimate of the local error for each cell in the mesh and are subsequently used to mark the mesh. In addition the global estimator can be used to determine when to stop iterating.

Mark

We have used two distinct marking strategies throughout the results section: the maximum strategy on the three-dimensional test cases and Dörfler strategy on the two-dimensional ones. We follow the presentation in [START_REF] Pfeiler | Dörfler marking with minimal cardinality is a linear complexity problem[END_REF]. In the maximum marking strategy [START_REF] Babuška | Feedback and adaptive finite element solution of onedimensional boundary value problems[END_REF], a cell is marked if its indicator is greater than a fixed fraction of the maximum indicator. More precisely, given a marking fraction θ ∈ (0, 1], the marked set M ⊂ T is the subset such that:

η bw,T ≥ θ max T ∈T η bw,T , ∀T ∈ T . (2.29)
In the Dörfler marking strategy [START_REF] Dörfler | A Convergent Adaptive Algorithm for Poisson's Equation[END_REF] (sometimes referred to as the equilibrated marking strategy) enough elements must be marked such that the sum of their estimators is larger than a fixed fraction of the total error. Given a marking fraction θ ∈ (0, 1], the marked set M is the subset with minimal cardinality #M such that

T ∈M η 2 bw,T ≥ θ T ∈T η 2 bw,T . (2.30)
We implement an O(N log N ) with N := #T complexity algorithm for finding the minimum cardinality set by sorting the indicators in decreasing order and finding the cutoff point such that eq. (2.30) is satisfied. Because of the ordering operation this set is guaranteed to have minimal cardinality. We note that recent work [START_REF] Hoare | Algorithm 65: find[END_REF][START_REF] Pfeiler | Dörfler marking with minimal cardinality is a linear complexity problem[END_REF] proposes a O(N ) complexity algorithm for finding the set with minimum cardinality.

Refine

We use two-dimensional and three-dimensional variants of the algorithm proposed in [START_REF] Plaza | Local refinement of simplicial grids based on the skeleton[END_REF], sometimes referred to as the Plaza algorithm. This algorithm works by subdividing the facets of each marked triangle or tetrahedron cell and then subdividing each triangle or tetrahedral cell so that it is compatible with the refinement on the facets. The algorithm has O(M ) complexity in the number of added mesh vertices M . This algorithm already exists in DOLFIN [START_REF] Logg | DOLFIN: Automated finite element computing[END_REF] and was used for the numerical results in [START_REF] Rognes | Automated Goal-Oriented Error Control I: Stationary Variational Problems[END_REF].

Goal-oriented adaptive mesh refinement

In many practical applications it is desirable to control the error in a specific quantity of interest, rather than the (global, i.e. across the entire domain Ω) energy norm [START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF]. In this section we show how the basic Bank-Weiser estimator can be used to control error in a goal functional, rather than in the natural norm. To do this, we use a weighted marking strategy proposed in [START_REF] Becker | Weighted Marking for Goal-oriented Adaptive Finite Element Methods[END_REF]. Let J : L 2 (Ω) → R be a given linear functional. Associated with J (u) and the primal problem eq. (2.2) is the dual or adjoint problem: Find the dual solution z ∈ H

1 D (Ω) such that (∇v, ∇z) = J (v), ∀v ∈ H 1 D (Ω).
(2.31)

The dual problem, like the primal problem, can also be approximated using the finite element method. Find z k ∈ V k such that

(∇v k , ∇z k ) = J (v k ) = (c, v k ) + (h, v k ) Γ , ∀v k ∈ V k . (2.32)
Using Galerkin orthogonality and Cauchy-Schwarz, it follows that

|J (u) -J (u k )| = |(∇(u -u k ), ∇z)| (2.33) = |(∇(u -u k ), ∇(z -z k ))| (2.34) ≤ ∇(u -u k ) ∇(z -z k ) , (2.35)
where the inequality holds due to Galerkin orthogonality.

Approximating the primal and dual errors ∇(u -u k ) and ∇(z -z k ) with any estimators η u and η z respectively, gives us an estimator for the error in the goal functional |J(u) -J(u k )| as the product of η u and η z , thanks to eq. (2.35):

η w := η u η z (2.36)
In addition, if η u and η z are reliable estimators i.e. if there exist two constants C u and C z only depending on the mesh regularity such that

∇(u -u k ) ≤ C u η u , and ∇(z -z k ) ≤ C z η z , (2.37)
then, η w is reliable as well

|J(u) -J(u k )| C u C z η w . (2.38)
Note that because the error in the goal functional is bounded by the product of two estimates, the element marking strategy must incorporate information from local indicators for both approximations to reduce the error on refinement. There are multiple strategies for doing this in the literature, see e.g. [START_REF] Mommer | A Goal-Oriented Adaptive Finite Element Method with Convergence Rates[END_REF]. We have chosen to implement the weighted goal-oriented (WGO) marking strategy from [START_REF] Becker | Weighted Marking for Goal-oriented Adaptive Finite Element Methods[END_REF]. The local WGO estimator is then defined as

η 2 w,T := η 2 z η 2 u + η 2 z η 2 u,T + η 2 u η 2 u + η 2 z η 2 z,T , ∀T ∈ T . (2.39)
The marking and refinement using η 2 w,T then follows in exactly the same manner as in the standard adaptive refinement strategy.

Extension to linear elasticity problems

Our implementation of the Bank-Weiser estimator can be directly applied to mixed formulations of (nearly-incompressible) linear elasticity problems using the results in [START_REF] Khan | Robust a posteriori error estimators for mixed approximation of nearly incompressible elasticity[END_REF]. In [START_REF] Ainsworth | A Posteriori Error Estimators for the Stokes and Oseen Equations[END_REF] a new a posteriori error estimator is introduced for mixed formulations of Stokes problems consisting in solving a local Poisson problem based on the local residuals on each cell. This estimator has been proved to be reliable and efficient in [START_REF] Ainsworth | A Posteriori Error Estimators for the Stokes and Oseen Equations[END_REF] under a saturation assumption. This assumption has been later removed in [START_REF] Liao | A simple yet effective a posteriori estimator for classical mixed approximation of Stokes equations[END_REF]. The reliability and efficiency of the estimator for mixed formulations of linear elasticity is proved in [START_REF] Khan | Robust a posteriori error estimators for mixed approximation of nearly incompressible elasticity[END_REF] without the need of a saturation assumption. In addition, they show that the estimator is robust in the incompressible limit.

Nearly-incompressible elasticity

We consider the problem of linear deformation of an isotropic elastic solid Ω using the Herrmann mixed formulation. We consider the stress tensor σ : Ω → R d×d , the strain tensor ε : Ω → R d×d , the load

f : Ω → R d which belongs to L 2 (Ω) d , the Dirichlet boundary data u D in H 1/2 (Γ D ) d ,
2µ (ε(u), ε(v)) -(p, div(v)) = (f , v) + (g, v) Γ N , ∀v ∈ H 1 D (Ω), (2.41a)
(q, div(u)) + 1 λ (p, q) = 0, ∀q ∈ L 2 (Ω). (2.41b)
The problem given by eqs. (2.41a) and (2.41b) admits a unique solution (see e.g. [START_REF] Khan | Robust a posteriori error estimators for mixed approximation of nearly incompressible elasticity[END_REF]). We introduce the finite element spaces

X D ⊂ H 1 D (Ω) and M ⊂ L 2 (Ω) such that X D := V 2 D d , (2.42)
and

M := V 1 . Let w be a discretization of u D ∈ X.
Considering the stable Taylor-Hood method of discretization, the mixed finite element approximation of eqs. (2.41a) and (2.41b) reads: find For a cell T and an edge E the residuals are defined by

u 2 ∈ X D with u 2 = w on Γ D and p 1 ∈ M such that 2µ (ε(u 2 ), ε(v 2 )) -(p 1 , div(v 2 )) = (f , v 2 ) + (g, v 2 ) , ∀v 2 ∈ X D , (2.43a)
(q 1 , div(u 2 )) + 1 λ (p 1 , q 1 ) = 0, ∀q 1 ∈ M. ( 2
R T := (f + div (2µε(u 2 )) -∇p 1 ) |T , (2.44a) r T := div(u 2 ) + 1 λ p 1 |T , (2.44b) R E =    1 2 (p 1 Id -2µε(u 2 )) n if E ∈ E I , 0 if E ∈ E D , g -(p 1 Id -2µε(u 2 ))n if E ∈ E N , (2.44c)
Here, once again we derive the a posteriori error estimator from these residuals and a local Poisson problem, following [START_REF] Khan | Robust a posteriori error estimators for mixed approximation of nearly incompressible elasticity[END_REF]. Let T be a cell of the mesh, the local Poisson problem read: find e T ∈ V bw T such that

2µ (∇e T , ∇v T ) T = (R T , v T ) T - E∈∂T (R E , v T ) E , ∀v T ∈ V bw T . (2.45)
The Poisson estimator is then defined by

η 2 p := T ∈T η 2 p,T , (2.46a) η 2 p,T := 2µ ∇e T 2 T + r T 2 T .
(2.46b) This estimator has been proved to be reliable and locally efficient in [START_REF] Khan | Robust a posteriori error estimators for mixed approximation of nearly incompressible elasticity[END_REF] as well as robust in the incompressible limit.

Results

We illustrate our implementation first on several two dimensional problems as Poisson problems with solutions of different regularities and with different boundary conditions. Then, we also look at examples of linear elasticity, and goal-oriented problems. We now treat a three dimensional example: a linear elasticity problem on a mesh inspired by a human femur bone. One can find another example of three dimensional application in [START_REF] Bulle | Removing the saturation assumption in Bank-Weiser error estimator analysis in dimension three[END_REF].

All the numerical results were produced within DOLFIN except the strong scaling tests in Section 2.9 which were performed using the DOLFINx version of our code.

We apply different adaptive refinement methods as presented in Section 2.5.1. For each method we perform the estimation step with a different estimator among the following: η res the residual estimator, defined in appendix A, η zz the Zienkiewicz-Zhu estimator, defined in appendix B. Note that we use the most basic version of the Zienkiewicz-Zhu estimator which is not defined for quadratic or cubic finite elements nor for linear elasticity problems, and consequently will be absent from the comparison in these cases (It is possible to extend the idea of the Zienkiewicz-Zhu estimator to higher-order polynomials via the definition of the Scott-Zhang interpolator, see [START_REF] Carstensen | Axioms of adaptivity[END_REF][START_REF] Scott | Finite element interpolation of nonsmooth functions satisfying boundary conditions[END_REF]). In addition we compare several versions of the Bank-Weiser estimator: the bubble Bank-Weiser estimator η b bw defined from the enriched bubble functions space and η k + ,k - bw for multiple choices of the fine and coarse spaces orders k + and k -.

For each one of the following test cases we will first give a comparison of all the refinement strategies by giving the efficiency of the a posteriori error estimator on the last mesh of the hierarchy, where the efficiency of an estimator η is defined as follows:

eff := η ε err , (2.47)
where ε err is a higher order approximation of the exact error computed either from the knowledge of the analytical solution or from a higher-order finite element method on a fine mesh.

Poisson problems

Returns -------

The error estimate on each cell of the mesh. """ mesh = u_h.function_space().mesh()

#

L-shaped domain

We consider a 2D L-shaped domain Ω = (-1, 1) 2 \ [-1, 0] 2 . We solve eq. (2.1) with f = 0, Γ D = Γ, u D given by the analytical solution defined below and Γ N = ∅. In polar coordinates, the exact solution is given by u exact (r, θ) = r 2/3 sin 2/3(θ + π/2) . The exact solution belongs to H 5/3-ε (Ω) for any ε > 0 and its gradient admits a singularity at the vertex of the reentrant corner [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF]Chapter 5]. L-shaped domains are widely used to test adaptive mesh refinement procedures [START_REF] Mitchell | A collection of 2D elliptic problems for testing adaptive grid refinement algorithms[END_REF]. In both linear and quadratic finite elements all the estimators reach an expected convergence rate (≈ -0.5 in the number of degrees of freedom for linear elements and ≈ -1 for quadratic elements). The choice of a posteriori error estimator is not critical for mesh refinement purposes, every estimator leading to a hierarchy of meshes on which the corresponding errors ε err are similar. For brevity we have not included the convergence plots of these results.

Linear elements. On fig. 2.2 we can see the initial mesh (top left) used to start the adaptive refinement strategies. Then, we can see the different refined meshes we obtain for k -> 1 largely underestimates it, leading to poor error approximations. Among the poor estimators, η 3,2 bw is surprisingly off for linear elements on this test case. This behavior seems to be specific to the L-shaped test cases with linear finite elements as we will see below.

Quadratic elements. As shown on fig. 2.4, the best estimator in terms of efficiency is η 2,0 bw which nearly perfectly matches the error ε err . We can also notice the very good efficiencies of η 4,2 bw and η 3,2 bw . Once again the Bank-Weiser estimators with k -> 2 drastically underestimate the error. We can notice that the residual estimator is less efficient as the finite element degree increases.

Mixed boundary conditions L-shaped domain

We solve eq. (2.1) on the same two-dimensional L-shaped boundary domain as in Section 2.7.2 but with different boundary conditions. We consider f = 0, Γ N = {(x, y) ∈ R 2 , x < 0, y = 0} and Γ D = Γ \ Γ N . The boundary data are given by g = 0 and u D = u exact = r 1/3 sin 1/3(θ + π/2) . The exact solution belongs to H 4/3-ε (Ω) for any ε > 0 and its gradient has a singularity located at the reentrant corner of Γ (see [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF]Chapter 5]). As before, each estimator is leading to a convergence rate close to the expected one (≈ -0.5 for linear elements, ≈ -1 for quadratic elements) and the choice of the estimator does not impact the quality of the mesh hierarchy.

Linear elements. First thing we can notice from fig. 2.5 is that the estimators efficiencies are quite different from those in fig. 2.3. Most of the Bank-Weiser estimator efficiencies have improved, except when k -> 1. The Zienkiewicz-Zhu estimator η zz is no longer the most efficient and has been outperformed by η 2,0 bw , η 2,1 bw and η 3,0 bw . The Bank-Weiser estimator η 3,2 bw still performs poorly as in fig. 2.3, while the residual estimator η res once again largely overestimates the error.

Quadratic elements. As for linear elements, the efficiencies in fig. 2 and other estimators on the last mesh of an adaptively refined hierarchy.

for the previous test cases, the Bank-Weiser estimators with k -> 2 are largely underestimating the error.

Boundary singularity

We solve eq. (2.1) on a two-dimensional unit square domain Ω = (0, 1) 2 with u = u exact on Γ D = Γ, (Γ N = ∅) and f chosen in order to have u(x, y) = u exact (x, y) = x α , with α 0.5.

In the following results we chose α = 0.7. The gradient of the exact solution u admits a singularity along the left boundary of Ω (for x = 0). The solution u belongs to H 6/5-ε for all ε > 0 [START_REF] Houston | Sobolev regularity estimation for hp-adaptive finite element methods[END_REF][START_REF] Mitchell | A collection of 2D elliptic problems for testing adaptive grid refinement algorithms[END_REF]. Consequently, the value of α determines the strength of the singularity and the regularity of u.

Due to the presence of the edge singularity, all the estimators are achieving a convergence rate close to -0.2 for linear elements. Moreover, this rate does not improve for higher-order elements (for brevity, the results for higher-order elements are not shown here). The low convergence rate shows how computationally challenging such a problem can be. Once again the choice of estimator is not critical for mesh refinement purposes.

Linear elements. The best estimator in terms of efficiency is η 2,1 bw which slightly over- and other estimators on the last mesh of an adaptively refined hierarchy.

k - k + 3 2
estimates the error, closely followed by η 4,2 bw underestimating the error as we can see on fig. 2.7. Unlike the previous test case, here the Zienkiewicz-Zhu estimator η zz grandly underestimates the error. The worst estimator is the residual estimator η res which gives no precise information about the error. We can notice that the poor performance of the estimator η 3,2 bw on the L-shaped test case does not reproduce here.

Quadratic elements. Again, fig. 2.8 shows that the best estimator is η 2,1 bw closely followed by η 2,0 bw and the bubble estimator η b bw . The residual estimator is getting worse as the finite element degree increases.

Goal-oriented adaptive refinement using linear elements

We solve the L-shaped domain problem as described in Section 2.7.2 but instead of controlling the error in the natural norm, we aim to control the error in the goal functional

J(u) = (c, u) with c a smooth bump function c(r) := ε -2 exp -1 r2 0 ≤ r2 < 1, 0 r2 ≥ 1. (2.48)
where r2 = ((x -x) /ε) 2 + ((y -ȳ) /ε) 2 , with ε ∈ R a parameter that controls the size of the bump function, and x ∈ R and ȳ ∈ R the position of the bumps function's center. We We use the goal-oriented adaptive mesh refinement methodology outlined in Section 2.5.2. We use a first-order polynomial finite element method for the primal and dual problem, and the Bank-Weiser error estimation procedure to calculate both η u and η z .

The 'exact' value of the functional J(u) was calculated on a very fine mesh using a fourth-order polynomial finite element space and was used to compute higher-order approximate errors for each refinement strategy.

The weighted goal-oriented strategy refines both the re-entrant corner and the broader region of interest defined by the goal functional. Relatively less refinement occurs in the regions far away from either of these important areas.

Figure 2.9 shows refined meshes after seven iterations of the weighted goal oriented method. We can see that the meshes are mainly refined in the re-entrant corner as well as in the region on the right top of it where the goal functional focuses. In fig. 2.10 we show the convergence curves of some of these adaptive strategies. For each strategy, η u = η z is the estimator specified in the legend and η w = η u η z . All the strategies we have tried led to very similar higher-order approximate errors. So for the sake of clarity we have replaced the approximate errors by an indicative line computed using a regression from the least squares method (lstsq error), leading to the line that fits the best the values of the different approximate errors. As we can see, these adaptive strategies are reaching an optimal convergence rate. Although it is also the case for all the other strategies we have tried, we do not show the other results for the sake of concision. In the left table of fig. 2.11 we show the efficiencies of the estimators η w where

η u = η z = η k + ,k - bw
. On the right table of fig. 2.11 we take η u = η z to be the estimators in the left column. As we can see on the efficiencies are not as good as in Section 2.7.2. The two best estimators are those derived from η 4,0 bw and η 4,2 bw . With η 3,0 bw , they are the only cases where the goal-oriented estimator η w is performing better than the goal-oriented estimator derived from the Zienkiewicz-Zhu estimator. The estimators η w derived from the bubble Bank-Weiser estimator as well as from the residual estimator are poorly overestimating the error. and other estimators on the last mesh of an adaptively refined hierarchy.

Linear elasticity problems

= (f 1 , f 2 ) with f 1 (x, y) = -2µπ 3 cos(πy) sin(πy) 2 cos(2πx) -1 , f 2 (x, y) = 2µπ 3 cos(πx) sin(πx) 2 cos(2πy) -1 . (2.49)
The corresponding exact solution of the linear elasticity problem reads u = (u 1 , u 2 ) with u 1 (x, y) = π cos(πy) sin 2 (πx) sin(πy), u 2 (x, y) = -π cos(πx) sin(πx) sin 2 (πy), (2.50)

the Herrmann pressure is zero everywhere on Ω. In each case we discretize this problem using the Taylor-Hood element and an initial Cartesian mesh and we apply our adaptive procedure driven by the Poisson estimator described in Section 2.5.3. We compare the Poisson estimators derived from different Bank-Weiser estimators and the residual estimator.

As before, all the refinement strategies are achieving an optimal convergence rate no matter the value of ν. fig. 2.12 shows the results for ν = 0.3. We notice that almost all the Poisson estimators derived from Bank-Weiser estimators have a very good efficiency. The best estimator in this case is η 2,0 bw closely followed by η 3,0 bw , η b bw and η 4,0 bw . Although the residual estimator still performs the worst, it is sharper than in all the previous test cases. As we can notice on fig. 2.13, all the estimators are robust with respect to the incompressibility constraint. All the efficiencies have slightly increased and some estimators (η 2,0 bw and η 3,0 bw ) that where a lower bound of the error previously are now an upper bound.

Human femur modeled using linear elasticity

In this test case we consider a linear elasticity problem on a domain inspired by a human femur bone 1 .

The goal of this test case is not to provide an accurate description of the behavior of the femur bone but to demonstrate the applicability of our implementation to 3D dimensional goal-oriented problem with large number of degrees of freedom: the linear elasticity problem to solve on the initial mesh, using Taylor-Hood element has 247,233 degrees of freedom while our last refinement step reaches 3,103,594 degrees of freedom. and other estimators on the last mesh of an adaptively refined hierarchy.

The 3D mesh for analysis is build from the surface model using the C++ library CGAL [START_REF] Alliez | 3D mesh generation[END_REF] via the Python front-end pygalmesh. The material parameters, namely the Young's modulus is set to 20 GPa and the Poisson's ratio to 0.42 (see e.g. [START_REF] Rupin | Experimental determination of Young modulus and Poisson ratio in cortical bone tissue using high resolution scanning acoustic microscopy and nanoindentation[END_REF]). In addition, the load is given by f = (0, 0, 0), the Dirichlet data by u D = 0 on Γ D Γ represented as the left dark gray region of the boundary in fig. 2.14 and g the traction data is defined as g = (0, 0, 0) on the center light gray region of the boundary and is constant on the right dark gray region of the boundary g = (-10 -7 , -10 -7 , 10 -6 ). The femur-shaped domain Ω as well as the initial and last meshes are shown in fig. 2.14. As we can see, the refinement occurs mainly in the central region of the femur, where the goal functional J focus. Some artifacts can be seen as stains of refinement in the central region due to the fact that we use the initial mesh as our geometry and on the left due to the discontinuity in the boundary conditions.

In fig. 2.15 the primal solution is given by the couple (u 2 , p 1 ) and the dual solution by (z 2 , κ 1 ). As we can notice and as expected, the weighted estimator η w converges twice as fast as the primal and dual estimators.

Strong scaling study

Finally, we provide results showing that our implementation scales strongly in parallel and that for a large-scale three-dimensional problem this error estimation takes significantly less time than the solution of the primal problem. In this section we use the new DOLFINx solver [START_REF] Habera | FEniCSX: A sustainable future for the FEniCS Project[END_REF] with the matching implementation of our algorithm.

We briefly discuss some aspects that are important for interpretation of the results. For a given cell the computation of the Bank-Weiser estimator requires geometry and solution data on the current cell and on all cells attached across its facets. So in a parallel computing context, cells located on the boundary of a partition require data from cells owned by another process. Both DOLFIN and DOLFINx support facet-mode ghosting where all data owned by cells on a partition boundary that share a facet are duplicated by the other process (ghost data). After the solution of the primal linear system the ghost data is updated between processes, which requires parallel communication. After this update, each process has a local copy of all of the data from the other rank needed to compute the Bank-Weiser estimator, and so the computation of the estimator is entirely local to a rank, i.e. without further parallel communication.

Because of this locality a proper implementation of this algorithm should demonstrate strong scaling performance. Furthermore, it would be desirable that the error estimation takes significantly less time than the solution of the primal problem even when using stateof-the art linear solution strategies. The results in this section demonstrate that this is indeed the case.

We solve eq. (2.1) where Ω is the unit cube [0, 1] 3 , Γ D = ∂Ω and Γ N = ∅. The data of this problem are given by f (x, y, z) = 12π 2 sin(2πx) sin(2πy) sin(2πz) and u D (x, y, z) = 0. These interior and boundary data are chosen such that the solution u of eq. (2.1) is given by u(x, y, z) = sin(2πx) sin(2πy) sin(2πz). We use continuous quadratic Lagrange finite elements and the Bank-Weiser error estimation is performed using the pair V 3 T /V 2 T . The primal linear system matrix and right-hand side vector are assembled using standard routines in DOLFINx. The resulting linear system is solved with PETSc [START_REF] Balay | PETSc Users Manual[END_REF] using the conjugate gradient method preconditioned with Hypre BoomerAMG algebraic multigrid [START_REF] Falgout | hypre: A Library of High Performance Preconditioners[END_REF].

The strong scaling study was carried out on the Aion cluster within the HPC facilities of the University of Luxembourg [START_REF] Varrette | Management of an academic HPC cluster: The UL experience[END_REF]. The Aion cluster is a Atos/Bull/AMD supercomputer composed of 318 compute nodes each containing two AMD Epyc ROME 7H12 processors with 64 cores per processor (128 cores per node). The nodes are connected through a Fast InfiniBand (IB) HDR 100Gbps interconnect in a 'fat-tree' topology. We invoke jobs using SLURM and ask for a contiguous allocation of nodes and exclusivity (no competing jobs) on each node. DOLFINx and PETSc are built using GCC 10.2.0 with Intel MPI and OpenBLAS. We use DOLFINx through its Python interface. The problem size is kept fixed at around 135 million degrees of freedom and the number of MPI ranks is increased from 128 (1 node, no interconnect communication) through to 2048 (16 nodes, interconnect communication) by doubling the number of nodes and ranks used in the previous computation.

In fig. 2.16 we show the results of the strong scaling study. We show wall time against MPI ranks and dof per rank for the primal linear system assembly, primal linear system solve, and the error estimation. For error estimation we are measuring steps 2 through 5 of Section 2.4.1. Both the solve and estimation scale almost perfectly down to around 65 thousand dof per rank. The primal system assembly does not scale as well as the estimation. This is because the primal system assembly is constrained by communication overheads and memory bandwidth, whereas the Bank-Weiser estimator computation is fully local and has much higher arithmetic intensity, so has not yet hit bandwidth limits of our system on the largest run. A further study (results not shown) using 96 MPI ranks per node yielded lower wall times and better strong scaling for primal linear system assembly, but the overall time for estimation and linear system solve increased and dominated any gains made in assembly. Comparing linear system assembly and solve with estimation time we can see that estimation is approximately one order of magnitude faster than solve time.

Conclusions

In this paper we have shown how the error estimator of Bank-Weiser, involving the solution of a local problem on a special finite element space, can be mathematically reformulated and implemented straightforwardly in a modern finite element software with the aid of automatic code generation techniques. Through a series of numerical results we have shown that the estimator is highly competitive in accurately predicting the global total error and in driving an adaptive mesh refinement strategy. Furthermore, the basic methodology and implementation for the Poisson problem can be extended to tackle more complex mixed discretizations of PDEs including nearly-incompressible elasticity or Stokes problems. We have also shown the (strong) scalability of our method when implemented in parallel and that the error estimation time is significantly lower than the primal solution time on a large problem.

Introduction

Fractional partial differential equations (FPDEs) are now applied in a wide range of fields [START_REF] Lischke | What is the fractional Laplacian? A comparative review with new results[END_REF] such as anomalous diffusion [START_REF] Bonito | Numerical methods for fractional diffusion[END_REF][START_REF] Carlson | A scalable framework for solving fractional diffusion equations[END_REF][START_REF] Defterli | Fractional Diffusion on Bounded Domains[END_REF][START_REF] Duo | A comparative study on nonlocal diffusion operators related to the fractional Laplacian[END_REF][START_REF] Nochetto | A PDE Approach to Fractional Diffusion in General Domains: A Priori Error Analysis[END_REF], electromagnetism and geophysical electromagnetism [START_REF] Bonito | Electroconvection of thin liquid crystals: Model reduction and numerical simulations[END_REF][START_REF] Weiss | Fractional Operators Applied to Geophysical Electromagnetics[END_REF], phase fluids [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF][START_REF] Akagi | Fractional Cahn-Hilliard, Allen-Cahn and porous medium equations[END_REF][START_REF] Du | Time-Fractional Allen-Cahn Equations: Analysis and Numerical Methods[END_REF], porous media [START_REF] Akagi | Fractional Cahn-Hilliard, Allen-Cahn and porous medium equations[END_REF][START_REF] Caputo | Models of flux in porous media with memory[END_REF], quasi-geostrophic flows [START_REF] Bonito | Numerical Simulations of Surface Quasi-Geostrophic Flows on Periodic Domains[END_REF] and spatial statistics [START_REF] Bolin | Numerical solution of fractional elliptic stochastic PDEs with spatial white noise[END_REF][START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF].

Fractional models can reproduce non-local behavior with a small number of parameters [START_REF] Atangana | Fractional Operators and Their Applications[END_REF][START_REF] Caputo | Linear Models of Dissipation whose Q is almost Frequency Independent-II[END_REF]. This non-locality, although useful from a modeling perspective, is a challenge for numerical methods since it naturally leads to large dense linear systems that can be computationally expensive to solve.

In the last decade various numerical methods have been derived in order to circumvent the main issues associated with the application of standard numerical methods to FPDEs, the two main ones being the non-locality leading to dense linear systems and, for some particular definitions of the fractional operator, the evaluation of singular integrals [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF][START_REF] Ainsworth | Towards an Efficient Finite Element Method for the Integral Fractional Laplacian on Polygonal Domains[END_REF].

We focus on discretization schemes based on finite element methods, other methods can be found e.g. in [START_REF] Aceto | Rational Krylov methods for functions of matrices with applications to fractional partial differential equations[END_REF][START_REF] Higham | An Improved Schur-Padé Algorithm for Fractional Powers of a Matrix and Their Fréchet Derivatives[END_REF][START_REF] Podlubny | Matrix approach to discrete fractional calculus II: Partial fractional differential equations[END_REF]. Among the methods addressing the above numerical issues, we can cite: methods to efficiently solve eigenvalue problems [START_REF] Carlson | A scalable framework for solving fractional diffusion equations[END_REF], multigrid methods for performing efficient dense matrix-vector products [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF][START_REF] Ainsworth | Towards an Efficient Finite Element Method for the Integral Fractional Laplacian on Polygonal Domains[END_REF], hybrid finite elementspectral schemes [START_REF] Ainsworth | Hybrid Finite Element-Spectral Method for the Fractional Laplacian: Approximation Theory and Efficient Solver[END_REF], Dirichlet-to-Neumann maps (such as the Caffarelli-Silvestre extension) [START_REF] Antil | A FEM for an Optimal Control Problem of Fractional Powers of Elliptic Operators[END_REF][START_REF] Caffarelli | An Extension Problem Related to the Fractional Laplacian[END_REF][START_REF] Chen | A PDE approach to fractional diffusion: A posteriori error analysis[END_REF][START_REF] Faustmann | Local convergence of the FEM for the integral fractional Laplacian[END_REF][START_REF] Nochetto | A PDE Approach to Fractional Diffusion in General Domains: A Priori Error Analysis[END_REF][START_REF] Stinga | Extension Problem and Harnack's Inequality for Some Fractional Operators[END_REF], semi-groups methods [START_REF] Cusimano | Discretizations of the Spectral Fractional Laplacian on General Domains with Dirichlet, Neumann, and Robin Boundary Conditions[END_REF][START_REF] Cusimano | Numerical approximations for fractional elliptic equations via the method of semigroups[END_REF][START_REF] Stinga | User's guide to the fractional Laplacian and the method of semigroups[END_REF], rational approximation methods [START_REF] Aceto | Rational Approximation to the Fractional Laplacian Operator in Reaction-Diffusion Problems[END_REF][START_REF] Harizanov | Numerical solution of fractional diffusion-reaction problems based on BURA[END_REF][START_REF] Hofreither | A unified view of some numerical methods for fractional diffusion[END_REF], Dunford-Taylor integrals [START_REF] Bolin | Numerical solution of fractional elliptic stochastic PDEs with spatial white noise[END_REF][START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF][START_REF] Bonito | The approximation of parabolic equations involving fractional powers of elliptic operators[END_REF][START_REF] Bonito | Numerical Approximation of Space-Time Fractional Parabolic Equations[END_REF][START_REF] Bonito | On sinc quadrature approximations of fractional powers of regularly accretive operators[END_REF][START_REF] Borthagaray | Linear and nonlinear fractional elliptic problems[END_REF][START_REF] Gavrilyuk | Data-sparse approximation to a class of operator-valued functions[END_REF][START_REF] Hofreither | A unified view of some numerical methods for fractional diffusion[END_REF] (which can be considered as particular examples of rational approximation methods) and reduced basis methods [START_REF] Danczul | A Reduced Basis Method For Fractional Diffusion Operators I[END_REF][START_REF] Danczul | A Reduced Basis Method For Fractional Diffusion Operators II. may[END_REF][START_REF] Dinh | Model reduction for fractional elliptic problems using Kato's formula[END_REF].

Although we focus exclusively on the spectral definition of the fractional Laplacian, there is no unique definition of the fractional power of the Laplacian operator. The three most frequently found definitions of the fractional Laplacian are: the integral fractional Laplacian, defined from the principal value of a singular integral over the whole space R d [12, 13, 67, 78, 109], the regional fractional Laplacian, defined by the same singular integral but over a bounded domain only [START_REF] Chen | The Dirichlet elliptic problem involving regional fractional Laplacian[END_REF][START_REF] Duo | A comparative study on nonlocal diffusion operators related to the fractional Laplacian[END_REF][START_REF] Fall | Regional fractional Laplacians: Boundary regularity[END_REF][START_REF] Mou | Interior Regularity for Regional Fractional Laplacian[END_REF] and the spectral fractional Laplacian, defined from the spectrum of the standard Laplacian over a bounded domain [START_REF] Ainsworth | Hybrid Finite Element-Spectral Method for the Fractional Laplacian: Approximation Theory and Efficient Solver[END_REF][START_REF] Antil | Fractional operators with inhomogeneous boundary conditions: analysis, control, and discretization[END_REF][START_REF] Banjai | Tensor FEM for Spectral Fractional Diffusion[END_REF][START_REF] Cusimano | Numerical approximations for fractional elliptic equations via the method of semigroups[END_REF][START_REF] Harizanov | Numerical solution of fractional diffusion-reaction problems based on BURA[END_REF][START_REF] Meidner | hp-Finite Elements for Fractional Diffusion[END_REF]. The different definitions are equivalent in the entire space R d , but this is no longer the case on a bounded domain [START_REF] Bonito | Numerical methods for fractional diffusion[END_REF][START_REF] Duo | A comparative study on nonlocal diffusion operators related to the fractional Laplacian[END_REF][START_REF] Kwaśnicki | Ten equivalent definitions of the fractional laplace operator[END_REF][START_REF] Lischke | What is the fractional Laplacian? A comparative review with new results[END_REF]. These definitions lead to significantly different mathematical problems associated with infinitesimal generators of different stochastic processes [START_REF] Duo | A comparative study on nonlocal diffusion operators related to the fractional Laplacian[END_REF][START_REF] Lischke | What is the fractional Laplacian? A comparative review with new results[END_REF].

Efficient methods for solving fractional problems typically rely on a combination of different discretization methods. For example, [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF], which is also the foundation of this work, combines a quadrature scheme for the Dunford-Taylor integral representation of the spectral fractional Laplacian with a standard finite element method in space. Both the quadrature scheme and the finite element method induce discretization errors. Each of these schemes is associated with its own discretization error. In order to achieve a solution to a given accuracy while avoiding wasted computational time, these errors need to be balanced.

A priori error estimation has been tackled for some definitions of the fractional Laplacian, such as the integral Laplacian [START_REF] Acosta | A Fractional Laplace Equation: Regularity of Solutions and Finite Element Approximations[END_REF][START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF][START_REF] Ainsworth | Hybrid Finite Element-Spectral Method for the Fractional Laplacian: Approximation Theory and Efficient Solver[END_REF][START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF][START_REF] Bonito | Numerical methods for fractional diffusion[END_REF][START_REF] Gimperlein | Space-time adaptive finite elements for nonlocal parabolic variational inequalities[END_REF] and the spectral fractional Laplacian [START_REF] Antil | A FEM for an Optimal Control Problem of Fractional Powers of Elliptic Operators[END_REF][START_REF] Antil | Fractional operators with inhomogeneous boundary conditions: analysis, control, and discretization[END_REF][START_REF] Banjai | Tensor FEM for Spectral Fractional Diffusion[END_REF][START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF][START_REF] Meidner | hp-Finite Elements for Fractional Diffusion[END_REF][START_REF] Nochetto | A PDE Approach to Fractional Diffusion in General Domains: A Priori Error Analysis[END_REF]. Unlike the standard Laplacian equation, solutions to the fractional Laplacian problems often exhibit strong boundary layers even for smooth data, particularly when the fractional power is low. These singularities lead to computational difficulties and have to be taken into account using, for example a priori geometric mesh refinement towards the boundary of the domain [START_REF] Acosta | A Fractional Laplace Equation: Regularity of Solutions and Finite Element Approximations[END_REF][START_REF] Banjai | Exponential convergence of hp FEM for spectral fractional diffusion in polygons[END_REF][START_REF] Bonito | Numerical approximation of the integral fractional Laplacian[END_REF][START_REF] Borthagaray | Linear and nonlinear fractional elliptic problems[END_REF][START_REF] Faustmann | Local convergence of the FEM for the integral fractional Laplacian[END_REF], or partition of unity enrichments [START_REF] Bordas | Partition of unity methods[END_REF]. We emphasize that [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] contains already an a priori error analysis in the L 2 norm for the combined rational sum finite element method that we use in this work.

A posteriori error estimation has also been considered in the literature on fractional equations. A simple residual based estimator is proposed for the integral fractional Laplacian in [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF]. A similar idea is used in the context of non-local variational inequalities in [START_REF] Gimperlein | Space-time adaptive finite elements for nonlocal parabolic variational inequalities[END_REF][START_REF] Nochetto | A posteriori error analysis for a class of integral equations and variational inequalities[END_REF]. Gradient-recovery based a posteriori error estimation has been developed in the context of fractional differential equations in [START_REF] Zhao | Adaptive finite element method for fractional differential equations using hierarchical matrices[END_REF]. In [START_REF] Bonito | Numerical methods for fractional diffusion[END_REF][START_REF] Chen | A PDE approach to fractional diffusion: A posteriori error analysis[END_REF] the authors present another estimator, based on the solution to local problems on cylindrical stars, for the integral fractional Laplacian discretized using the Caffarelli-Silvestre extension. A weighted residual estimator is derived in [START_REF] Faustmann | Quasi-optimal convergence rate for an adaptive method for the integral fractional laplacian[END_REF] in the same context.

To our knowledge, no a posteriori error estimation method has been derived for the spectral fractional Laplacian, discretized using the rational approximation approach of [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF].

Contribution

The main contribution of this work is the derivation of a novel a posteriori error estimator for the combined rational finite element approximation of the spectral fractional Laplacian. It is a natural a posteriori counterpart to the a priori results developed in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF].

Our work starts with the quadrature rule for the Dunford-Taylor integral proposed in the seminal work [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF]. This method, and other rational approximation-based discretization methods, decompose the original fractional problem into a set of independent parametric non-fractional problems. From this point we develop an associated set of independent non-fractional a posteriori error estimation problems. We compute the Bank-Weiser hierarchical estimators [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF] of the error between each non-fractional parametric problem solution and its finite element discretization, then the fractional problem discretization error is estimated by the sum of the parametric contributions via the rational approximation.

Our method leads to a fully local and parallelizable solution technique for the spectral fractional Laplacian with computable L 2 error. Our method is valid for any finite element degree (however, for the sake of brevity we do not show results with higher degree finite elements) and for one, two and three dimensional problems [START_REF] Bulle | Hierarchical a posteriori error estimation of Bank-Weiser type in the FEniCS Project[END_REF].

We implement our method in DOLFINx [START_REF] Habera | FEniCSX: A sustainable future for the FEniCS Project[END_REF], the new problem solving environment of the FEniCS Project [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF]. A simple demonstration implementation is included in the supplementary material. We show numerical results demonstrating that the estimator can correctly reproduce the a priori convergence rates derived in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF].

Our newly developed error estimator is then used to steer an adaptive mesh refinement algorithm, resulting in improved convergence rates for small fractional powers and strong boundary layers.

Motivation

Given a fractional power s in (0, 1) and a rational approximation Q κ s (λ) of the function λ -s , it is possible to construct a semi-discrete approximation u κ of the solution u to a fractional Laplace equation as a weighted sum of solutions (u l ) l to non-fractional parametric problems. Then, a fully discrete approximation of u is obtained by discretizing the parametric solutions (u l ) l using a finite element method.

An a posteriori error estimator is then computed as the weighted sum of the Bank-Weiser estimators of the error between each u l and its finite element discretization. As we will see in the following, the resulting numerical scheme is simple and its implementation in code is straightforward. Furthermore it maintains the appealing embarrassingly parallel nature of rational approximation schemes [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF][START_REF] Gavrilyuk | Data-sparse approximation to a class of operator-valued functions[END_REF][START_REF] Harizanov | Numerical solution of fractional diffusion-reaction problems based on BURA[END_REF].

We remark on why we have chosen to use the Bank-Weiser type error estimator, as opposed to one of the many other error estimation strategies, e.g. explicit residual, equilibrated fluxes, or recovery-type estimators (see [START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF][START_REF] Carstensen | Estimator Competition for poisson Problems[END_REF] and references therein). In the case of fractional powers of the Laplacian operator, the resulting set of parametric problems consists of singularly-perturbed reaction-diffusion equations. It has been proven in [START_REF] Verfürth | Robust a posteriori error estimators for a singularly perturbed reactiondiffusion equation[END_REF] that the Bank-Weiser estimator is robust with respect to the coefficients appearing in these parametric problems when the error is measured in the natural norm. To our knowledge, no such robustness has been established for the L 2 -norm for the Bank-Weiser estimator. Nevertheless, our numerical experiments indicate that this does appear to be the case. Moreover, the Bank-Weiser estimator can be straightforwardly applied to higher-order finite element methods and higher-dimension problems. In addition, its computational stencil is highly local which is particularly appealing for three-dimensional problems see e.g. [START_REF] Bulle | Hierarchical a posteriori error estimation of Bank-Weiser type in the FEniCS Project[END_REF].

In this work we focus on error estimation in the L 2 norm, the estimation of the error in the 'natural' fractional norm is the topic of ongoing work. For simplicity, we only consider fractional powers of the Laplacian with homogeneous Dirichlet boundary conditions.

Problem statement

For any subset ω of Ω we denote L 2 (ω) the space of square integrable functions on ω and (•, •) ω its usual inner product. Let H 1 (ω) be the Sobolev space of functions with first order weak derivatives in L 2 (ω). The space H 1 (ω) is endowed with the usual inner product

(∇•, ∇•) L 2 (ω) + (•, •) L 2 (ω)
. We will omit the dependence in ω in the subscripts when ω = Ω. We will make use of the notation ∂v/∂n := ∇v • n for the normal derivative of a smooth enough function v. We denote H 1 0 (Ω) the subspace of functions in H 1 (Ω) with a zero trace on Γ.

We consider the family of eigenfunctions {ψ i } ∞ i=1 ⊂ H 1 0 (Ω) of the standard Laplacian operator with uniform zero Dirichlet boundary condition on Ω as well as the corresponding family of eigenvalues {λ i } ∞ i=1 . We assume the Laplacian eigenvalues are sorted in increasing order and we assume λ 0 ∈ R is a lower bound of the spectrum

λ 0 λ 1 • • • λ i λ i+1 • • • (3.1)
The family {ψ i } ∞ i=1 is an orthonormal basis of L 2 (Ω). For s in (0, 1) we introduce the spectral fractional Sobolev space H s and its natural norm

H s := v ∈ L 2 (Ω), ∞ i=1 λ s i (v, ψ i ) 2 < ∞ , v 2 H s := ∞ i=1 λ s i (v, ψ i ) 2 . (3.2)
Especially, for 0 s 1 we have

H 1 0 (Ω) = H 1 (Ω) ⊆ H s (Ω) ⊆ L 2 (Ω) =: H 0 (Ω)
and the norm • H s coincide with • L 2 when s = 0 and with |•| H 1 when s = 1.

The spectral fractional Laplacian

Let s be a real number in (0, 1) and f be a given function in L 2 (Ω). We consider the following fractional Laplacian problem: we look for a function u such that

(-∆) s u = f in Ω, u = 0 on Γ. (3.3)
The solution u of eq. (3.3) is defined using the spectrum of the standard Laplacian [START_REF] Antil | Fractional operators with inhomogeneous boundary conditions: analysis, control, and discretization[END_REF] u :=

∞ i=1 λ -s i (f, ψ i ) ψ i . (3.4)
For a data f in L 2 (Ω), the solution u belongs to H 2s [START_REF] Antil | Fractional operators with inhomogeneous boundary conditions: analysis, control, and discretization[END_REF] and

u H 2s = f L 2 . (3.5)
Using the spectrum of the Laplacian, we can derive a weaker formulation of eq. (3.3). Multiplying eq. (3.3) by test functions v in H s (Ω) and integrating over Ω gives

Ω (-∆) s uv = Ω f v, ∀v ∈ H s (Ω). (3.6)
Now, considering the expansions of u and v in the basis {ψ i } ∞ i=1 we have

((-∆) s u, v) = ∞ i=1 λ s i u i v i = ∞ i,j=1 λ s/2 i u i λ s/2 j v j (ψ i , ψ j ) = (-∆) s/2 u, (-∆) s/2 v .
Then, the solution u to eq. (3.6) satisfies

(-∆) s/2 u, (-∆) s/2 v = (f, v) , ∀v ∈ H s (Ω). (3.7)
The formulation eq. (3.7) is weaker since it is valid for u in H s (Ω) only while for eq. (3.3) the solution u belongs to H 2s (Ω). If we assume the solution to eq. (3.7) to be in H 2s then the formulations eq. (3.7) and eq. (3.3) are equivalent.

Rational approximation

Our method relies on rational approximations of the real function λ → λ -s for s in (0, 1) and λ λ 0 for some fixed λ 0 > 0. We are particularly interested in an example provided in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF]. This example is based on the following expression [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] λ -s = 2 sin(πs) π

+∞

-∞ e 2sy 1 + e 2y λ -1 dy.

(3.8)

Then, the rational approximation is obtained from eq. (3.8) by discretizing the integral on the right-hand side with a trapezoidal quadrature rule,

λ -s Q κ s (λ) := 2 sin(πs) πκ N (κ) l=-M (κ) e 2slκ 1 + e 2lκ λ -1 , (3.9)
where κ > 0 is the fineness parameter and M (κ) := π 2 4sκ 2 , and N (κ) :=

π 2 4(1 -s)κ 2 , (3.10)
where • is the ceiling function.

This particular scheme has some advantages compared to other rational methods. The coefficients (e 2slκ ) N l=-M and (e 2lκ ) N l=-M are very easy to compute in comparison with methods based on e.g. best uniform rational approximations (BURA) (see [START_REF] Aceto | Rational Approximation to the Fractional Laplacian Operator in Reaction-Diffusion Problems[END_REF][START_REF] Harizanov | Numerical solution of fractional diffusion-reaction problems based on BURA[END_REF][START_REF] Hofreither | A unified view of some numerical methods for fractional diffusion[END_REF][START_REF] Hofreither | An algorithm for best rational approximation based on barycentric rational interpolation[END_REF]). This scheme is also among the most efficient as shown in recent comparison studies (see [START_REF] Hofreither | A unified view of some numerical methods for fractional diffusion[END_REF][START_REF] Vabishchevich | Approximation of a fractional power of an elliptic operator[END_REF]). Various other examples of rational approximations can be found e.g. in [START_REF] Aceto | Rational Approximation to the Fractional Laplacian Operator in Reaction-Diffusion Problems[END_REF][START_REF] Aceto | Rational approximations to fractional powers of self-adjoint positive operators[END_REF][START_REF] Gavrilyuk | Data-sparse approximation to the operator-valued functions of elliptic operator[END_REF][START_REF] Harizanov | Numerical solution of fractional diffusion-reaction problems based on BURA[END_REF][START_REF] Vabishchevich | Approximation of a fractional power of an elliptic operator[END_REF]. We want to highlight again that the error estimation scheme developed later can be derived in the same manner regardless of the choice of the rational approximation, as long as it leads to a set of well-posed non-fractional parametric problems.

It has been shown in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] that Q κ s converges uniformly to λ -s at an exponential rate as κ → 0. Especially, the approximation error is bounded by

λ -s -Q κ s (λ) ε s (κ), ∀λ λ 0 , ∀κ > 0, (3.11)
with

ε s (κ) = 2 sin(πs) π 1 2s + 1 2(1 -s)λ 0 1 1 -e -π 2 /(2κ) + 1 e -π 2 /(2κ) .
(3.12) Asymptotically, ε s (κ) behaves like e -π 2 /(2κ) as κ → 0.

Discretization

We combine the rational approximation eq. (3.9) with a finite element method to derive a fully discrete approximation of the solution u to eq. (3.7).

Rational semi-discrete approximation

From eq. (3.9) we can derive semi-discrete approximations of the solution u to eq. (3.7) by considering

u κ := 2 sin(πs) πκ N (κ) l=-M (κ) e 2slκ u l , (3.13)
where the functions {u l } N l=1 are solutions to the parametric problems: for each l in -M, N , find u l in H 1 0 such that

(u l , w) + e 2slκ (∇u l , ∇w) = (f, w) ∀w ∈ H 1 0 . (3.14)
It has been proved in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] that the semi-discrete approximation u κ converges to u in L 2 (Ω) at the same speed as Q κ s (λ) converges to λ -s . More precisely,

u -u κ L 2 ε s (κ) f L 2 , ∀κ > 0. (3.15)
where ε s (κ) is defined in eq. (3.12). We can deduce from eq. (3.15) the following two important points. Firstly, the rational approximation u κ converges to u exponentially fast in κ. Therefore, it does not constitute a bottleneck in the rate of convergence when combined with a finite element method to obtain a fully discrete approximation. Secondly, the right-hand side of eq. (3.15) is technically an a posteriori estimation of the rational discretization error since ε s (κ) and f L 2 can be calculated almost entirely using a priori known data. The only parameter that is not so easily computable in ε s is λ 0 , a lower bound of the spectrum of the Laplacian on Ω. The bound ε s (κ) can be optimized by taking λ 0 = λ 1 but given its exponential convergence rate, ε s (κ) will not drastically deteriorates if λ 0 < λ 1 . Moreover, precise guaranteed lower bounds for λ 1 could be obtained following e.g. [START_REF] Cancès | Guaranteed and Robust a Posteriori Bounds for Laplace Eigenvalues and Eigenvectors: Conforming Approximations[END_REF][START_REF] Carstensen | Axioms of adaptivity[END_REF].

Finite element discretization

In order to get a fully discrete approximation of u, we use a finite element method to discretize the parametric problems eq. (3.14). Although it is not mandatory, we use the same mesh and same finite element space for all the parametric problems. We discuss this choice, and possible alternative strategies, in section 3.7.

Let T be a mesh on the domain Ω, composed of cells T = {T }, facets E = {E} (we call facets the edges in dimension two and the faces in dimension three), and vertices. The mesh T is supposed to be regular, in Ciarlet's sense: h T /ρ T γ, ∀T ∈ T , where h T is the diameter of a cell T , ρ T the diameter of its inscribed ball, and γ is a positive constant fixed once and for all. The subset of facets that are not coincident with the boundary Γ (called interior facets) is denoted E I . Let n + and n -in R d be the outward unit normals to a given edge as seen by two cells T + and T -incident to a common edge E. The space of polynomials of order p on a cell T is denoted P p (T ) and the continuous Lagrange finite element space of order p on the mesh T is defined by

V p := v p ∈ H 1 (Ω), v p|T ∈ P p (T ) ∀T ∈ T . (3.16)
We denote V p 0 the finite element space composed by functions of V p vanishing on the boundary Γ. For a given index l, the finite element discretization of eq. (3.14) reads: for each l in -M, N , find u l,p in V p 0 such that (u l,p , v p ) + e 2lκ (∇u l,p , ∇v p ) = (f, v p ) , ∀v p ∈ V p 0 .

(3.17)

Then, combining eq. (3.13) with eq. (3.17) we can give a fully discrete approximation of the solution to eq. (3.7)

u ≈ u κ,p := 2 sin(πs) πκ

N (κ) l=-M (κ)
e 2slκ u l,p .

(3.18)

The computation of u κ,p is summarized in the top part of fig. 3.1.

Finite element discretization error analysis

According to what we have seen in section 3.4.2, the rational approximation error, characterized by u -u κ L 2 converges exponentially fast. Consequently, we will consider this error to be negligible and assume that the rational scheme Q κ s is precise enough (i.e. κ is small enough) so that u u κ .

Our goal is to bound the discretization error in the

L 2 norm u -u κ,p L 2 u κ -u κ,p L 2 . (3.20)
Since for any s ∈ (0, 1), the discrepancy u -u κ,p belongs to H s (Ω) ⊂ L 2 (Ω), the error can be measured in the L 2 norm for any value of the fractional power s.

Heuristics

Let us start with some heuristics motivating the derivation of our a posteriori error estimator. The main idea is to derive a function e bw κ,T that locally represents the discretization error in the solution to the fractional problem (u κ -u κ,p ) |T on a cell T of the mesh. Thanks to the rational approximation we notice that

(u κ -u κ,p ) |T = 2 sin(πs) πκ N l=-M e 2slκ (u l -u l,p ) |T . (3.21)
So we can use the framework proposed by Bank and Weiser in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF] to derive solutions e bw l,T such that e bw l,T

(u l -u l,p ) |T , ∀l ∈ -M, N , ∀T ∈ T . (3.22)
We obtain e bw κ,T using the rational approximation sum We would like to emphasize that the Bank-Weiser estimator is not the only possible choice. In fact, the Bank-Weiser estimator could be replaced with another estimator based on the solves of local problems, such as e.g. the one used in [START_REF] Nochetto | A PDE Approach to Fractional Diffusion in General Domains: A Priori Error Analysis[END_REF].

A posteriori error estimation

Let us now derive our a posteriori error estimation method more precisely. As mentioned in the last subsection, this estimator is based on a hierarchical estimator computed from the solves of local Neumann problems on the cells and introduced for the first time by Bank and Weiser in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF].

Let T be a cell of the mesh. We make use of the following local finite element spaces

V p T := v p,T ∈ P p (T ), v p,T = 0 in (Ω \ T ) ∪ (T ∩ ∂Ω) . (3.25)
Let us now consider two non-negative integers p + and p -such that p + > p - 0 and

L T : V p + T -→ V p - T
the local Lagrange interpolation operator. We introduce the local Bank-Weiser space, defined by

V bw T := ker(L T ) = v p + ,T ∈ V p + T , L T (v p + ,T ) = 0 . (3.26)
The local parametric Bank-Weiser problem associated to the parametric problems eq. (3.14) and eq. (3.17) reads

T e bw l,T v bw T + e 2slκ T ∇e bw l,T • ∇v bw T = T r l,T v bw T + 1 2 E∈∂T E J l,E v bw T , ∀v bw T ∈ V bw T (3.27)
where r l,T and J l,T are defined as follow:

r l,T := f |T -u l,p |T + e 2lκ ∆u l,p |T , and J l,T := e 2lκ ∂u l,p ∂n E .

(3.28)

The solution e bw l,T in V bw T is the local parametric Bank-Weiser solution. More details about the computation and implementation of the Bank-Weiser solutions can be found in [START_REF] Bank | Some A Posteriori Error Estimators for Elliptic Partial Differential Equations[END_REF][START_REF] Bulle | Hierarchical a posteriori error estimation of Bank-Weiser type in the FEniCS Project[END_REF].

Then, we derive the local fractional Bank-Weiser solution by summing the local parametric Bank-Weiser solutions into the rational approximation sum The local fractional Bank-Weiser estimator is then defined as the L 2 norm of this local solution η bw κ,T := e bw κ,T L 2 (T ) .

(3.30) The global fractional Bank-Weiser estimator is then defined by

u -M u -M +1 u 0 u N -1 u N Rational sum u κ u -M,p u -M +1,p u 0,p u N -1,p u N,p u κ,p
η bw κ 2 := T ∈T η bw κ,T 2 .
(3.31)

Adaptive refinement

One of the main applications of a posteriori error estimation is to drive adaptive mesh refinement algorithms. When the error is unevenly spread across the mesh, refining uniformly is a waste of computational resources leading to sub-optimal convergence rates in the number of degrees of freedom. This problem is compounded for computationally expensive problems like fractional problems. Moreover, it is known that fractional problems often show a boundary layer behavior, the discretization error is consequently large in a localized region near the boundary [START_REF] Acosta | A short FE implementation for a 2d homogeneous Dirichlet problem of a fractional Laplacian[END_REF][START_REF] Borthagaray | Local Energy Estimates for the Fractional Laplacian[END_REF][START_REF] Vabishchevich | Approximation of a fractional power of an elliptic operator[END_REF]. This problem has been tackled using graded meshes that are refined near the boundary based on a priori or a posteriori considerations [START_REF] Bonito | Numerical methods for fractional diffusion[END_REF][START_REF] Chen | A PDE approach to fractional diffusion: A posteriori error analysis[END_REF][START_REF] Gimperlein | Space-time adaptive finite elements for nonlocal parabolic variational inequalities[END_REF][START_REF] Meidner | hp-Finite Elements for Fractional Diffusion[END_REF]. As expected, the use of graded meshes improves the convergence of the methods. Adaptive refinement algorithms are based on the loop

• • • -→ Solve -→ Estimate -→ Mark -→ Refine -→ • • •
In this work we are concerned with developments in the modules solve and estimate. We are using totally standard approaches, namely the Dörfler algorithm [START_REF] Dörfler | A Convergent Adaptive Algorithm for Poisson's Equation[END_REF] for the mark module and the Plaza-Carey algorithm [START_REF] Plaza | Local refinement of simplicial grids based on the skeleton[END_REF] for the refine module. Rational approximation methods have the advantage of being fully parallelizable due to the independence of the parametric problems from each other. Similarly, the local a posteriori error estimation method we have presented earlier is also parallelizable since the Frac. power 0.1 0.3 0.5 0.7 0.9 Num. param. prob. 408 176 149 176 408 Table 3.1: Number of parametric problems solved for each fractional power. parametric problem at each iteration. More importantly, this would mean summing parametric finite element solutions coming from different and possibly non-nested meshes. Properly addressing this question is beyond the scope of this study. Nonetheless, we give some hints the numerical section section 3.9.1.

Implementation

We have implemented our method using the DOLFINX finite element solver of the FEniCS Project [START_REF] Alnaes | The FEniCS Project Version 1.5[END_REF]. Each parametric subproblem is submitted to a batch job queue. A distinct MPI communicator is used for each job. We use a standard first-order Lagrange finite element method and the resulting linear system is solved using the conjugate gradient method preconditioned using BoomerAMG from HYPRE [START_REF] Falgout | hypre: A Library of High Performance Preconditioners[END_REF] via the interface in PETSc [START_REF] Balay | PETSc Users Manual[END_REF]. To compute the Bank-Weiser error estimator for each subproblem we use the methodology outlined in [START_REF] Bulle | Hierarchical a posteriori error estimation of Bank-Weiser type in the FEniCS Project[END_REF] and implemented in the FEniCSx-EE package [START_REF] Bulle | FEniCS Error Estimation (FEniCS-EE)[END_REF]. For every subproblem the computed solution and error estimate is written to disk in HDF5 format. A final step, running on a single MPI communicator, reads the solutions and error estimates for all subproblems, computes the quadrature sums using axpy operations, defines the marked set of cells to be refined using the Dörfler algorithm [START_REF] Dörfler | A Convergent Adaptive Algorithm for Poisson's Equation[END_REF], and finally refines the mesh using the Plaza-Carey algorithm [START_REF] Plaza | Local refinement of simplicial grids based on the skeleton[END_REF].

A more complex implementation using a single MPI communicator split into several sub-communicators would remove the necessity of reading and writing the solution and error estimate for each subproblem to and from disk. However, in practice the cost of computing the parametric solutions massively dominates all other costs.

Numerical results

First, we need to choose the value of κ in order to guarantee that the rational approximation error is negligible. From eq. (3.15), we know a bound that depends on s, λ 0 and f L 2 . However, in all our test cases we know that λ 0 = λ 1 = 1 is a lower bound for the spectrum of the Laplacian and the data f is always chosen such that f L 2 = 1. It turns out that taking κ = 0.26 ensures that u -u κ L 2 10 -8 , no matter the choice of s ∈ (0, 1). This choice leads to a different number of parametric problems to solve for each fractional power, these numbers are detailed in table 3.1.

When analytical solutions are known, we provide the efficiency indices of the Bank-Weiser estimator, defined by η bw κ / u κ -u κ,1 L 2 .

Two-dimensional product of sines test case

We solve eq. ( 3.3) on the square Ω = (0, π) 2 with data f (x, y) = (2/π) sin(x) sin(y). The analytical solution to this problem is given by u(x, y) = 2 -s (2/π) sin(x) sin(y). Moreover, the analytical solutions to the parametric problems eq. (3.14) are also known u l (x, y) = (1+ 2 e 2lκ ) -1 (2/π) sin(x) sin(y). The problem is solved on a hierarchy of structured (triangular) meshes. For this test case the solution u shows no boundary layer behavior, therefore adaptive refinement cannot improve the convergence rate. Consequently we only perform uniform refinement on this case. As we can see on fig. 3.3, the Bank-Weiser estimator tends to be very accurate when the mesh is fine enough when s = 0.3 and s = 0.7. In fact, its accuracy is robust with respect to the fractional power. The efficiency indices are computed by taking the average of the ratios for the five last meshes of the hierarchy and are shown for various fractional powers in table 3.3. Theorem 4.3 from [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] gives a convergence rate for the finite element scheme depending on the elliptic regularity index α of the Laplacian over Ω, on the fractional power and on the regularity index δ of the data f . Since Ω is convex the elliptic "pick-up" regularity index α can be taken to be 1 [START_REF] Bonito | Numerical methods for fractional diffusion[END_REF] and since f is infinitely smooth the coefficient δ can be taken as large as wanted. Consequently, Theorem 4.3 in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] predicts a convergence rate of dof -1 for this test case. The convergence rates we measure in practice, shown in table 3.2, are coherent with this prediction. These rates are computed from a linear regression fit on the values obtained on the five last meshes of the hierarchy.

Parametric problems discretization error

Since we know the analytical solutions to the parametric problems in this case, it is possible to compute the exact parametric discretization errors u l -u l,1 L 2 , for each l ∈ -M, N . It is possible then to investigate the consequences of using the same mesh for all the parametric problems. In fig. 3.4 we have plotted the exact parametric errors after five steps of (uniform) refinement. As we can notice, the same mesh leads to a wide range of parametric errors values, especially for fractional powers s close to 1. These errors are particularly low for high values of the index l, when the diffusion part of the operator is dominant. However, when l becomes less than zero, i.e. when the reaction part is dominant the mesh seems to have an equal effect on the parametric errors. As expected these results suggest that the method can be optimized by using different meshes depending on l. In particular, coarser meshes would be sufficient for high values of l. These results are obtained for uniform refinement, further investigations deserve to be carried out for adaptive refinement.

As we explained earlier, using a different hierarchy of meshes for each parametric problem may be computationally advantageous, at the expense of ease of implementation. Several hierarchies of meshes would need to be stored and, in the case of adaptive mesh refinement, interpolation between possibly non-nested meshes would be required in order to compute the fractional solution u. To avoid these complications when adaptive refinement is used, we propose the following:

1. use the same hierarchy of meshes for all the parametric problems but not the same Frac. power 0.1 0.3 0.5 0.7 0.9 Estimator -0.92 -0.97 -0.99 -1.00 -1.00 Exact error -1.00 -1.00 -1.00 -1.00 -0.94 Table 3.2: Two-dimensional product of sines test case: convergence rates of the Bank-Weiser estimator and of the exact error for various fractional powers. mesh. Some parametric problems might be solved on coarser meshes from the hierarchy and others on finer ones. This would allow to keep only one hierarchy of meshes stored in memory. Moreover, it would avoid the interpolation between non-nested meshes, since meshes from the same hierarchy are always nested.

2. selectively refine the mesh hierarchy: estimate the error globally for each parametric problem (this can be done using the local parametric Bank-Weiser solutions) and mark the parametric problems for which a finer mesh is required, using e.g. a marking algorithm similar to Defiler's marking strategy.

Frac. power 0.1 0.3 0.5 0.7 0.9 Est. eff. index 0.86 1.16 1.08 0.96 0.78 Frac. power 0.1 0.3 0.5 0.7 0.9 Estimator -0.56 -0.60 -0.63 -0.65 -0.66 Exact error -0.69 -0.69 -0.69 -0.69 -0.69 Table 3.4: Three-dimensional product of sines test case: convergence rates of the Bank-Weiser estimator and of the exact error for various fractional powers.

Three-dimensional product of sines test case

This test case is the three-dimensional equivalent of the last test case. We solve eq. (3.3) on the cube Ω = (0, π) 3 with data f (x, y, z) = (2/π) 3/2 sin(x) sin(y) sin(z). The analytical solution to this problem is given by u(x, y, z) = 3 -s (2/π) 3/2 sin(x) sin(y) sin(z). The problem is solved on a hierarchy of uniformly refined Cartesian (tetrahedral) meshes. As for the two-dimensional case, the solution u shows no boundary layer behavior and adaptive refinement is not required. For the same reasons as for the two-dimensional case, Theorem 4.3 from [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] predicts a convergence rate of dof -2/3 for the finite element scheme. fig. 3.5 shows the values of the Bank-Weiser estimator and of the exact error (computed from the knowledge of the analytical solution) for s = 0.3 and s = 0.7. As in the twodimensional case, the efficiency indices are relatively robust with respect to the fractional powers. They are shown for various fractional powers in table 3.5 and are computed by taking the average of the indices from the three last meshes of the hierarchy. As we can see, the Bank-Weiser estimator efficiency indices for this three-dimensional case are not as good as in the two-dimensional case. We have already observed this behavior for nonfractional problems [START_REF] Bulle | Removing the saturation assumption in Bank-Weiser error estimator analysis in dimension three[END_REF]. We can notice that the convergence rates, given in table 3.4, are coherent with the predictions of Theorem 4.3 from [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF]. The convergence rates are computed from a linear regression on the values computed from the three last meshes of the hierarchy. Frac. power 0.1 0.3 0.5 0.7 0.9 Est. eff. index 2.12 3.20 3.08 2.77 2.45 Table 3.5: Three-dimensional product of sines test case: efficiency indices of the Bank-Weiser estimator for various fractional powers.

Two-dimensional checkerboard test case

We solve the problem introduced in the numerical results of [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF]. We consider a unit square Ω = (0, 1) 2 with data f : Ω → R given for all (x 1 , x 2 ) ∈ Ω by

f (x 1 , x 2 ) = 1, if (x 1 -0.5)(x 2 -0.5) > 0, 0, otherwise. (3.32)
The data f ∈ H 1/2-ε (Ω) for all ε > 0. So in Theorem 4.3 of [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] the index δ < 1/2 and since Ω is convex again α can be chosen equal to 1. Then, the predicted convergence rate (for uniform refinement) is ln(

√ dof)dof -β with β = 1, if s > 3 4 , s + 1 4 , otherwise. (3.33)
The predicted (if we omit the logarithmic term) and calculated convergence rates for different choices of s are given in table 3.6. As we can see on this table, the convergence rates for the Bank-Weiser estimator is globally coherent with the predictions. fig. 3.6 shows that adaptive refinement improves the convergence rate for small fractional powers. This is expected, the deterioration in the convergence rate is due to the boundary layer behavior of the solution that is getting stronger as the fractional power decreases. When the fractional κ when uniform refinement is performed (dashed light lines) and when adaptive refinement is performed (solid lines).

power is close to 1, the solution behaves like the solution to a non-fractional problem for which adaptive refinement is no longer needed. This can be seen on fig. 3.7, after 10 steps of adaptive refinement, the mesh associated to fractional power s = 0.9 is almost uniformly refined while the meshes associated to s = 0.5 and s = 0.1 show strongly localized refinement. This explains why in fig. 3.6 we barely see any improvement in the convergence rate when the mesh is adaptively refined compared to uniformly refined when s 0.7.

Three-dimensional checkerboard test case

This test case is the three-dimensional version of the above checkerboard problem. We solve eq. (3.3) on the unit cube Ω = (0, 1) 3 , with data f such that

f (x 1 , x 2 , x 3 ) =      1, if (x 1 -0.5)(x 2 -0.5) > 0 and (x 3 -0.5) < 0, 1, if (x 1 -0.5)(x 2 -0.5) < 0 and (x 3 -0.5) > 0, -1, otherwise. (3.34)
The finite element solution u 1 and the corresponding mesh after six steps of adaptive refinement are shown in fig. 3.8 for the fractional power s = 0.5. As for the two-dimensional case, f ∈ H 1/2-ε (Ω) for all ε > 0. Consequently, once again Theorem 4. 3.6: Two-dimensional checkerboard test case: convergence slopes of the Bank-Weiser estimator for uniform refinement and for adaptive refinement compared to the values predicted by [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] for various fractional powers. 3.7: Three-dimensional checkerboard test case: convergence slopes of the Bank-Weiser estimator for uniform refinement and for adaptive refinement compared to the values predicted by [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] for various fractional powers. a convergence rate (for uniform refinement) equal to ln dof 1/3 dof -2β/3 with β given by eq. (3.33).

Once again, if we omit the logarithmic term, the predicted and calculated convergence rates are given in table 3.7. As in the two-dimensional case, the convergence rates of the Bank-Weiser estimator are globally coherent with the predictions and the boundary layer behavior becomes stronger as the fractional power decreases leading to poorer convergence rates. section 3.9.4 shows the values of the Bank-Weiser estimator for uniform and adaptive refinement and for several fractional powers. κ when uniform refinement is performed (dashed light lines) and when adaptive refinement is performed (solid lines).

Concluding remarks

In this work we presented a novel a posteriori error estimation method for the spectral fractional Laplacian. This method benefits from the embarrassingly parallel character of both the Bank-Weiser error estimator and the rational approximation methods, thus keeping the appealing computational aspects of the underlying methodology in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF]. Here are two important points we want to make to conclude this paper. First, the Bank-Weiser estimator seems to be equivalent to the L 2 exact error at least when structured meshes are used and when the solution u is smooth. Second, adaptive refinement methods drastically improves the convergence rate compared to uniform refinement for fractional powers close to 0.

Finally, we give some future directions that we think are worth considering. More numerical tests could be performed, especially for higher order elements and/or using variants of the Bank-Weiser error estimator as considered in [START_REF] Bulle | Hierarchical a posteriori error estimation of Bank-Weiser type in the FEniCS Project[END_REF].

We would like also to study the derivation of an algorithm that allows to use different meshes to discretize the parametric problems in order to save computational time, as explained in section 3.9.1.

The a posteriori error estimation of the error in the "natural" norm of the problem i.e. the spectral fractional norm defined in eq. (3.2) is another extension of this work that is worth to consider.

The replacement of the Bank-Weiser estimator by an anisotropic a posteriori error estimator would improve the convergence rate even further in case of boundary layers, see e.g. [START_REF] Banjai | Tensor FEM for Spectral Fractional Diffusion[END_REF][START_REF] Faustmann | Local convergence of the FEM for the integral fractional Laplacian[END_REF],

Another interesting extension would be to test our method on fractional powers of other kinds of elliptic operators, following [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF], on another definition of the fractional Laplacian operator [START_REF] Bonito | Numerical approximation of the integral fractional Laplacian[END_REF] and/or other boundary conditions, following [START_REF] Antil | Fractional operators with inhomogeneous boundary conditions: analysis, control, and discretization[END_REF].

Chapter 4

Convergence of a rational approximation method

This chapter is original and not based on any submitted or published articles.

Abstract

We propose a novel convergence result for a particular rational approximation method based on a Dunford-Taylor integral introduced in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF]. Thanks to this new result we improve the convergence result [66, Theorem 3.2] of a semi-discretization method of the solution to a fractional Laplacian equation based on this rational approximation. We show the exponential convergence rate of the method while alleviating the data regularity assumption. Finally, we propose a conjecture on the existence of a rational approximation method allowing to further improve this convergence result and provide a numerical study suggesting that the method from [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] does not satisfy the assumptions of the conjecture.

Introduction

In [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF], a method to solve fractional Laplacian equations is derived, based on a Dunford-Taylor integral representation of the function λ → λ s (for s ∈ (0, 1)). The discretization consists in the combination of a quadrature rule, leading to a rational approximation, with a finite element method. The fractional problem is then reformulated into a family of nonfractional parametric reaction-diffusion problems that can be discretized and solved independently using the finite element method.

Originally proposed for fractional power of elliptic operators, this method has been extended to other kind of equations such as, parabolic equations involving fractional elliptic operators [START_REF] Bonito | The approximation of parabolic equations involving fractional powers of elliptic operators[END_REF], space-time fractional parabolic equations [START_REF] Bonito | Numerical Approximation of Space-Time Fractional Parabolic Equations[END_REF], equations based on regularly accretive operators [START_REF] Bonito | On sinc quadrature approximations of fractional powers of regularly accretive operators[END_REF] and integral fractional Laplacian equations [START_REF] Bonito | Numerical approximation of the integral fractional Laplacian[END_REF].

A convergence result for the scalar rational approximation method is established in [60, Lemma 3.4], showing the exponential convergence rate of the method. From this lemma, convergence results for the semi-discrete rational approximation of the solution to a fractional Laplacian equation are proposed in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF]Theorem 3.5] and in [START_REF] Bonito | Numerical approximation of fractional powers of regularly accretive operators[END_REF]Theorem 7.1] for the L 2 norm and for a fractional Sobolev norm respectively. These results hold under the assumption that the data of the fractional equation is regular enough to belong to a fractional order Sobolev space.

Contribution

In this chapter we propose a novel convergence result (Lemma 1) for the scalar rational approximation method introduced in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF]. From this lemma a novel convergence result for the rational approximation of the solution u in a fractional Sobolev norm is derived (Theorem 2) . This result shows the exponential convergence rate of the method introduced in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] under a weaker regularity assumption on the input data, compared to [61, Theorem 7.1]. Finally, we make the conjecture of the existence of a rational approximation method satisfying a stronger convergence assumption whose would improve our results and simplify the proofs.

Problem statement

Let Ω be an bounded open domain of R d (d = 1, 2 or 3), with polygonal/polyhedral boundary denoted by Γ := ∂Ω. We denote L 2 (Ω) the Lebesgue space of square integrable functions on Ω and (•, •) and • L 2 its usual inner product and norm. We denote H 1 (Ω) the usual Sobolev space of order 1 and H 1 0 (Ω) the subspace of functions with a trace vanishing on Γ. We denote

(∇•, ∇•) and |•| H 1 the usual inner product and norm of H 1 0 (Ω). Let {(ψ i ) ∞ i=1 , (λ i ) ∞ i=1 } ⊂ L 2 (Ω) × R +, *
be the spectrum of the standard Laplacian operator with uniform zero Dirichlet boundary condition on Ω, where the eigenvalues (λ i ) ∞ i=1 are sorted in non-decreasing order. The family {ψ i } ∞ i=1 is an Hilbert (orthonormal) basis of L 2 (Ω). For s in [-1, 1] we introduce the spectral fractional Sobolev space [START_REF] Bonito | Numerical methods for fractional diffusion[END_REF]

H s (Ω) := v = ∞ i=1 v i ψ i , v i ∈ R ∀i ∈ N * : ∞ i=1 λ s i v 2 i < ∞ , (4.1)
of natural norm

v 2 H s := ∞ i=1 λ s i v 2 i . (4.2)
In particular, when s ∈ [0, 1] we have the following spaces inclusion

H s (Ω) ⊆ L 2 (Ω) and the coefficients v i of v ∈ H s (Ω) are given by v i = (v, ψ i ) for i ∈ N * . For s ∈ [0, 1], the space H -s (Ω) is the dual of H s (Ω) [65]. When s ∈ [-1, 1
] the space H s (Ω) can be characterized in terms of other fractional Sobolev spaces, see e.g. [START_REF] Antil | Fractional operators with inhomogeneous boundary conditions: analysis, control, and discretization[END_REF][START_REF] Bonito | Numerical methods for fractional diffusion[END_REF][START_REF] Lischke | What is the fractional Laplacian? A comparative review with new results[END_REF]. In particular, when -1 s 1 we have the following sequence of spaces inclusions

H 1 (Ω) ⊆ H s (Ω) ⊆ H -1 (Ω). (4.3)
In addition, for s = 0 (resp. s = 1),

H 0 (Ω) = L 2 (Ω) (resp. H 1 (Ω) = H 1 0 (Ω)) and the norm • H s coincide with the norm • L 2 (resp. |•| H 1 ).

The spectral fractional Laplacian

Let s be a real number in (0, 1) and f be a given function in

H δ (Ω) for δ ∈ [-1, 1], of coefficients (f i ) +∞ i=1 ⊂ R such that f = ∞ i=1 f i ψ i .
We consider the following fractional Laplacian problem: we seek a function u such that

(-∆) s u = f in Ω, u = 0 on Γ. (4.4)
The solution u of eq. (4.4) is defined using the spectrum of the standard Laplacian [START_REF] Antil | Fractional operators with inhomogeneous boundary conditions: analysis, control, and discretization[END_REF],

u := ∞ i=1 λ -s i f i ψ i . (4.5)
Then u ∈ H 2s+δ (Ω) and in particular,

u H 2s+δ = f H δ . (4.6)
We can derive an equivalent formulation of eq. (4.4). Multiplying eq. ( 4.4) by test functions v in H s (Ω) and integrating over Ω gives

Ω (-∆) s uv = Ω f v, ∀v ∈ H s (Ω). (4.7)
Now from the expansions of u and v in the basis {ψ i } ∞ i=1 we obtain

((-∆) s u, v) = ∞ i=1 λ s i u i v i = ∞ i,j=1 λ s/2 i u i λ s/2 j v j (ψ i , ψ j ) = (-∆) s/2 u, (-∆) s/2 v . (4.8)
Then, the solution u to eq. (4.7) satisfies

(-∆) s/2 u, (-∆) s/2 v = (f, v) , ∀v ∈ H s (Ω). (4.9)
Conversely, if u satisfy eq. (4.9), then using the decompositions of u, v and f in the basis (ψ i ) +∞ i=1 we obtain

u i = λ -s i f i , ∀i = 1, • • • , +∞, (4.10)
which is equivalent (by definition) to eq. (4.4).

Rational semi-discrete approximation

The method to derive a semi-discrete approximation to the solution u to eq. (4.9) relies on a rational approximation to the real function λ → λ -s for s in (0, 1) and λ λ 0 for some fixed λ 0 > 0. In [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF], the following integral representation is derived from the Balakrishnan's formula [START_REF] Balakrishnan | Fractional powers of closed operators and the semigroups generated by them[END_REF] λ -s = 2 sin(πs) π

+∞

-∞ e 2sy 1 + e 2y λ -1 dy.

When this integral is discretized using a rectangular quadrature rule, it is possible to derive the following rational approximation

λ -s Q N s (λ) := 2κ sin(πs) π N l=-N e 2slκ 1 + e 2lκ λ -1 , (4.12)
where κ = 1/ √ N . Various other examples of rational approximations can be found e.g. in [START_REF] Aceto | Rational Approximation to the Fractional Laplacian Operator in Reaction-Diffusion Problems[END_REF][START_REF] Aceto | Rational approximations to fractional powers of self-adjoint positive operators[END_REF][START_REF] Gavrilyuk | Data-sparse approximation to the operator-valued functions of elliptic operator[END_REF][START_REF] Harizanov | Numerical solution of fractional diffusion-reaction problems based on BURA[END_REF][START_REF] Vabishchevich | Approximation of a fractional power of an elliptic operator[END_REF]. Comparisons between different rational approximation schemes can be found in [START_REF] Čiegis | Parallel solvers for fractional power diffusion problems[END_REF][START_REF] Čiegis | A Comparison of Accuracy and Efficiency of Parallel Solvers for Fractional Power Diffusion Problems[END_REF][START_REF] Hofreither | A unified view of some numerical methods for fractional diffusion[END_REF][START_REF] Vabishchevich | Approximation of a fractional power of an elliptic operator[END_REF].

From eq. ( 4.12) we can derive the semi-discrete approximation to the solution u to eq. (4.9). We denote u Q N s the semi-discrete approximation to u, defined by

u Q N s := 2κ sin(πs) π N l=-N e 2slκ u l , (4.13)
where the functions {u l } N l=-N are the solutions to the parametric problems: for each l in -N, N , find u l in H 1 0 such that

(u l , w) + e 2lκ (∇u l , ∇w) = (-∆) s/2 u, (-∆) s/2 w = (f, w) , ∀w ∈ H 1 0 . (4.14)
The function u Q N s is not fully discrete. A fully discrete approximation to u could be obtained e.g. by discretizing the parametric problems eq. (4.14) using finite element methods.

Convergence of the rational approximation

In this section we study the convergence of the rational approximation defined in eq. (4.12). In particular, we show a convergence result slightly different from [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF]Lemma 3.4] where we explicit a dependence of the upper bound in λ. From this novel convergence result we deduce a new convergence result of the semi-discrete approximation u Q N s defined in eq. (4.13) to the solution u. This convergence result extends the result [64, Theorem 7.1].

Convergence of the scalar rational approximation

The result in [60, Lemma 3.4] shows that Q N s (λ) converges uniformly to λ -s at an exponential rate. In this study we establish a different convergence property for the rational approximation, given by the following lemma.

Lemma 1. Let s be in (0, 1), σ in [0, s) and λ 0 > 0. Then,

λ -s -Q N s (λ) λ -σ ε σ (N ), ∀λ ∈ [λ 0 , +∞), (4.15)
where ε σ (N ) is given by

ε σ (N ) := 2 sin(πs) π   π √ 2λ σ-s 0 e -(π 2 /2) √ N sin(πs) 1 -e -(π 2 /2) √ N + λ σ-1 0 e 2(s-1) √ N 2(1 -s) + e 2(σ-s) √ N 2(s -σ)   . (4.16) Especially, ε σ (N ) = O N →+∞ e 2 max(s-1,σ-s) √ N . (4.17)
Proof. We proceed as in [60, Section 3.3] and base our proof on Theorems 2.20 and 2.21 of [START_REF] Lund | Sinc methods for quadrature and differential equations[END_REF]. Let s ∈ (0, 1) and λ 0 > 0. Let D be the infinite strip of the complex plane defined by

D := {z ∈ C, |Im(z)| < π/4}. (4.18)
For λ λ 0 and z ∈ D, we denote g λ (z) := e 2sz 1 + e 2z λ -1 . According to eq. (4.11) and eq. (4.12), the rational approximation error can be written

λ -s -Q N s (λ) = 2 sin(πs) π +∞ -∞ g λ (y) dy -κ N l=-N g λ (lκ) , (4.19)
where κ := 1/ √ N . Thus, we have

λ -s -Q N s (λ) 2 sin(πs) π +∞ -∞ g λ (y) dy -κ +∞ l=-∞ g λ (lκ) (4.20)
+ κ +∞ l=N +1 g λ (lκ) + κ -N -1 l=-∞ g λ (lκ) .
Let us now bound the three terms on the right hand side. We denote

B 1 := +∞ -∞ g λ (y) dy -κ +∞ l=-∞ g λ (lκ) , (4.21)
B 2 := κ +∞ l=N +1 g λ (lκ) , (4.22)
B 3 := κ -N -1 l=-∞ g λ (lκ) . (4.23)
We start with B 1 . By [184, Theorem 2.20], if g λ is analytic in D and, for some 0 a < 1, satisfy

π/4 -π/4 |g λ (x + iη)|dη = O x→±∞ (|x| a ), (4.24)
and

N (g λ , D) := lim η→(π/4) - +∞ -∞ |g λ (x -iη)| + |g λ (x + iη)| dx < +∞, (4.25)
then we have

B 1 = +∞ -∞ g λ (y) dy -κ +∞ l=-∞ g λ (lκ) N (g λ , D) 2 sinh(π 2 /(4κ)) e -π 2 /(4κ) . (4.26)
Let us check that the function g λ satisfies all the hypotheses. First, for any λ > 0 the function g λ is analytic in D as a product of analytic functions in D. Second, let us check eq. (4.24) for x > 0. The function g λ can be bounded as follow

|g λ (x + iη)| = e 2sx |1 + e 2(x+iη) λ| -1 e 2(s-1)x λ -1 0 . (4.27)
Thus, since s < 1, for any 0 a < 1 we have

π/4 -π/4 |g λ (x + iη)|dη π/4 -π/4 e 2(s-1)x λ -1 0 dη = π 2 e 2(s-1)x λ -1 0 = O x→+∞ (|x| a ). (4.28)
Equation (4.24) for x < 0 is satisfied if we notice that

|g λ (x + iη)| e 2sx . (4.29)
Then, since s > 0, for any 0 a < 1 we have

π/4 -π/4 |g λ (x + iη)|dη π/4 -π/4 e 2sx dη = π 2 e 2sx = O x→-∞ (|x| a ). (4.30)
Finally, let us verify eq. (4.25). To do so, we compute this limit. The integrand in eq. (4.25) is given by

|g λ (x -iη)| + |g λ (x + iη)| = e 2s(x-iη) 1 + e 2(x-iη) λ -1 (4.31)
+ e 2s(x+iη) 1 + e 2(x+iη) λ -1

.

On the one hand we have, e 2s(x-iη) = e 2sx = e 2s(x-iπ/4) and e 2s(x+iη) = e 2sx = e 2s(x+iπ/4) . (

And on the other hand,

1 + e 2(x-iη) λ 2 = Re 1 + e 2(x-iη) λ 2 + Im 1 + e 2(x-iη) λ 2 = 1 + e 2x cos(-2η)λ 2 + e 2x sin(-2η)λ 2 = 1 + 2 e 2x cos(-2η)λ + e 2x λ 2 cos(-2η) 2 + sin(-2η) 2 = 2 e 2x cos(-2η)λ + 1 + e 4x λ 2 (4.33)
So for |η| π/4 we have

1 + e 2(x-iη) λ 1 + e 2(x-iπ/4) λ , (4.34)
and similarly

1 + e 2(x+iη) λ 1 + e 2(x+iπ/4) λ . and since λ > 0 we have

1 + e 2(x±iπ/4) λ 2 = 2 e 2x cos(±π/2) + 1 + e 4x λ 2 = 1 + e 4x λ 2 1 2 1 + e 2x λ 2 . So, +∞ -∞ |g λ (x -iπ/4)| + |g λ (x + iπ/4)| dx = +∞ -∞ e 2sx 1 + e 2(x-iπ/4) λ -1 + 1 + e 2(x+iπ/4) λ -1 dx 2 √ 2 +∞ -∞ e 2sx 1 + e 2x λ -1 dx. (4.37)
Using eq. (4.11) we get

+∞ -∞ |g λ (x -iπ/4)| + |g λ (x + iπ/4)| dx 2 √ 2 π 2 sin(πs) λ -s (4.38)
Since the integrand of eq. (4.25) is bounded by eq. (4.36) uniformly in η and since this bound is integrable eq. (4.38), by the dominated convergence theorem we have

N (g λ , D) = +∞ -∞ |g λ (x -iπ/4)| + |g λ (x + iπ/4)| dy. (4.39)
Then, if we take σ ∈ [0, s] by eq. (4.38) and since λ λ 0 we obtain

N (g λ , D) 2 √ 2 π 2 sin(πs) λ -s π √ 2 sin(πs) λ σ-s 0 λ -σ . (4.40)
Consequently, by eq. (4.26) for σ ∈ [0, s] we have a bound for the first of the three terms of eq. (4.20)

B 1 πλ σ-s 0 √ 2 sin(πs) sinh(π 2 /(4κ)) λ -σ e -π 2 /(4κ) . (4.41)
Using the definition of the hyperbolic sine we have,

e -π 2 /(4κ) sinh(π 2 /(4κ)) = 2 e -π 2 /(4κ) e π 2 /(4κ) -e -π 2 /(4κ) = 2 e -π 2 /(2κ) 1 -e -π 2 /(2κ) , (4.42)
so,

B 1 λ -σ π √ 2λ σ-s 0 e -π 2 /(2κ) sin(πs) 1 -e -π 2 /(2κ) (4.43)
We now bound the second term B 2 . We first establish the following useful inequality for ν in [0, 1], y ∈ R and λ > 0,

1 + e 2y λ -1 λ -ν e -2νy . (4.44)
The inequality eq. (4.44) is equivalent to

1 + e 2y λ -λ ν e 2νy 0 ∀y ∈ R. (4.45)
However, if we take x = λ e 2y then if x 1 we have x x ν so

1 + x -x ν 0. (4.46)
Now if 1 > x then 1 > x ν and 1 + x -x ν > 0. λ -ν e 2(s-ν)y . (

If we use eq. (4.48) with ν = 1 and the fact that λ 0 λ we have for any s in (0, 1) and any σ in [0, s],

B 2 = κ +∞ l=N +1 |g λ (lκ)| λ -1 κ +∞ l=N +1 e 2(s-1)lκ λ -σ λ σ-1 0 κ e 2(s-1)(N +1)κ +∞ l=0 e 2(s-1)lκ . (4.49)
Since s < 1, the geometric series on the right hand side is convergent and its sum is

+∞ l=0 e 2(s-1)lκ = 1 1 -e 2(s-1)κ . (4.50)
If we recall κ = 1 √ N and using the fact that e x

1-e x 1 -x for all x = 0, we have

B 2 λ -σ λ σ-1 0 κ e 2(s-1)N κ e 2(s-1)κ 1 -e 2(s-1)κ λ -σ λ σ-1 0 e 2(s-1)N κ κ 2(1 -s)κ λ -σ λ σ-1 0 2(1 -s) e 2(s-1)/κ . (4.51)
Bound of the last term B 3 . For s in (0, 1), using eq. (4.48) with ν = σ in [0, s) we have

B 3 = κ -N -1 l=-∞ |g λ (lκ)| = κ +∞ l=N +1 |g λ (-lκ)| = κ +∞ l=N +1 e -2slκ 1 + e -2lκ λ -1 λ -σ κ +∞ l=N +1 e 2(σ-s)lκ . (4.52)
Once again, the geometric series on the right hand side is convergent, of sum

+∞ l=N +1 e 2(σ-s)lκ = 1 1 -e 2(σ-s)κ . (4.53)
We can bound B 3 the same way we did for eq. (4.51) to obtain

B 3 λ -σ 1 2(s -σ) e 2(σ-s)/κ . (4.54)
Then, using eqs. (4.43), (4.51) and (4.54) in eq. (4.20) we obtain

λ -s -Q N s (λ) λ -σ ε σ (N ), (4.55)
with,

ε σ (N ) := 2 sin(πs) π π √ 2λ σ-s 0 e -π 2 /(2κ) sin(πs) 1 -e -π 2 /(2κ) + λ σ-1 0 e 2(s-1) √ N 2(1 -s) + e 2(σ-s) √ N 2(s -σ) . (4.56)
In particular,

ε σ (N ) = O N →+∞ e -2 min{1-s,s-σ} √ N . (4.57)
When σ = 0 we recover a result that is very similar to [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF]Lemma 3.4]. In [60, Remark 3.1] a rational approximation sum with an improved convergence rate is derived from Q N s by considering a non-symmetric sum. A similar thing can be achieved here.

In the proof of Lemma 1 if we replace Q N s by the following non-symmetric rational approximation sum depending on a fineness parameter 0 < κ 1,

Q κ s (λ) := 2κ sin(πs) π M (κ) l=-M (κ) e 2slκ 1 + e 2lκ λ -1 , (4.58)
with

M (κ) = π 2 4(1 -s)κ 2 ,
and

M (κ) = π 2 4(s -σ)κ 2 , (4.59)
it is possible to obtain the following upper bound

λ -s -Q κ s (λ) λ -σ ε σ (κ), ∀λ ∈ [λ 0 , +∞), (4.60)
where,

εσ (κ) = 2 sin(πs) π π √ 2λ σ-s 0 sin(πs) 1 -e -(π 2 /2)/κ + λ σ-1 0 2(1 -s) + 1 2(s -σ) e -(π 2 /2)/κ , (4.61)
In particular,

εσ (κ) = O κ→0 e -(π 2 /2)/κ . (4.62)
In comparison, if we explicit the dependence of ε σ in κ = 1/ √ N in eq. (4.57), we have

ε σ (κ) = O κ→0 e -2 min{1-s,s-σ}/κ . (4.63)
Note that due to the definition of M in eq. (4.59), the number of terms in the sum eq. (4.58) gets larger and larger as σ tends to s. Moreover, εσ and M become infinite when σ = s.

Rational approximation error analysis

The norm • H s is a natural choice to measure the approximation error. In the following theorem, we give a bound for the error between a function v and its semi-discrete approximation given by the rational sum Q N s defined in eq. (4.12).

Theorem 2. Let s ∈ (0, 1), σ ∈ [0, s] and v be a function of H 3s-2σ (Ω). Let in addition {v l } N l=-N be the solutions to the following problems: for each l ∈ -N, N , find

v l ∈ H 1 0 (Ω) such that (v l , w) + e 2lκ (∇v l , ∇w) = (-∆) s/2 v, (-∆) s/2 w , ∀w ∈ H 1 0 (Ω), (4.64)
and

v Q N s := 2κ sin(πs) π N l=-N e 2slκ v l , (4.65)
where κ = 1/ √ N . Then, lim

N →+∞ v -v Q N s H s = 0. (4.66)
Moreover, if σ < s we have

v -v Q N s H s ε σ (N ) v H 3s-2σ , (4.67)
where ε σ (N ) is defined in Lemma 1.

Proof. The proof follows the plan here below:

1. we first reduce this problem to a scalar function approximation problem and we use Lemma 1 for the case σ < s, 2. we then treat the remaining case σ = s independently.

1. Expanding the functions v, {v l } N l=-N and v Q N s in the basis of eigenfunctions of the Laplacian we have

v = +∞ i=1 v i ψ i , v l = +∞ i=1 v l,i ψ i , and v Q N s = +∞ i=1 v Q N s ,i ψ i . (4.68)
From eq. (4.64) we deduce for each l in -N, N and each i in 1, +∞

v l,i = 1 + e 2lκ λ i -1 λ s i v i , (4.69)
thus,

v Q N s ,i = Q N s (λ i )λ s i v i . (4.70)
Then, using eqs. (4.69) and (4.70) we have

v -v Q N s 2 H s = +∞ i=1 λ s i (v i -v Q N s ,i ) 2 (4.71) +∞ i=1 1 -Q N s (λ i )λ s i 2 λ s i v 2 i (4.72) +∞ i=1 λ -s i -Q N s (λ i ) 2 λ 3s i v 2 i . (4.73)
Now, using Lemma 1 we have, for any σ ∈ [0, s)

v -v Q N s 2 H s ε σ (N ) 2 +∞ i=1 λ 3s-2σ i v 2 i ε σ (N ) 2 v 2 H 3s-2σ . (4.74)
Consequently, when v ∈ H 3s-2σ (Ω),

lim N →+∞ v -v Q N s H s = 0, (4.75)
at an exponential rate.

2. The bound of Lemma 1 is not valid when σ = s because of the bound on the term B 3 , defined in eq. (4.23). Thus, we have to find another bound for B 3 .

First of all, for y 0 and λ λ 0 we have g λ (-y)

then if N M B 3 = κ -N -1 l=-∞ |g λ (lκ)| = κ +∞ l=N +1 |g λ (-lκ)| = κ M +1 l=N +1 |g λ (-lκ)| + κ +∞ l=M +2 |g λ (-lκ)| = S 1 + S 2 . (4.79)
The sum S 1 can be bounded using the same arguments as for eq. (4.51)

S 1 κ M +1 l=N +1 e 2(1-s)lκ λ κλ -1 e 2(1-s)(N +1)κ 1 -e 2(1-s)(M -N +1)κ 1 -e 2(1-s)κ κλ -1 e 2(1-s)κ 1 -e 2(1-s)κ e 2(1-s)N κ -e 2(1-s)(M +1)κ λ -1 κ 2(s -1)κ e 2(1-s)N κ -e 2(1-s)(M +1)κ . Now, since M 1 2κ ln(λ) then e 2(1-s)M κ λ 1-s . So S 1 λ -1 1 2(s -1) e 2(1-s)M κ e 2(1-s)κ e 2(1-s)(N -M -1)κ -1 λ -s 1 2(1 -s) e 2(1-s)κ 1 -e 2(1-s)(N -M -1)κ λ -s 1 2(1 -s) e 2(1-s)κ . (4.80)
For S 2 , using the arguments of eq. (4.51) we obtain

S 2 κ +∞ l=M +2
e -2slκ e -2s(M +1)κ 1 2s .

However, since M is the largest integer smaller than 1 2κ ln(λ), we have M + 1 > 1 2κ ln(λ). Then,

S 2 e -2s(M +1)κ 1 2s λ -s 1 2s . ( 4
B 3 λ -s 1 2(1-s) e 2(1-s)/ √ N + 1 2s if N 1 2 ln(λ) 2 λ -s 1 2s e -2sk(λ)/ √ N if N > 1 2 ln(λ) 2 , (4.85)
where k(λ) = N -1 2 ln(λ) 2 .

In summary, if σ = s and if

N 1 2 ln(λ) 2 then λ -s -Q N s (λ) 2 sin(πs) π (B 1 + B 2 + B 3 ) λ -s 2 sin(πs) π   π √ 2 e -(π 2 /2) √ N sin(πs) 1 -e -(π 2 /2) √ N + λ s-1 0 2(1 -s) e 2(s-1)
√ N (4.86)

+ 1 2(1 -s) e 2(1-s)/ √ N + 1 2s . If N > 1 2 ln(λ) 2 then λ -s -Q N s (λ) λ -s 2 sin(πs) π   π √ 2 e -(π 2 /2) √ N sin(πs) 1 -e -(π 2 /2) √ N + λ s-1 0 2(1 -s) e 2(s-1)
√ N (4.87)

+ 1 2s e -2sk(λ)/ √ N .
Now, if we introduce the index I min (N ) defined as

I min (N ) := min i ∈ 1, +∞ , λ i exp 2 √ N , (4.88)
we have i < I min (N ) ⇐⇒ N > 1 2 ln(λ i ) 2 . Then,

+∞ i=1 λ -s i -Q N s (λ i ) 2 λ 3s i v 2 i = I min (N )-1 i=1 λ -s i -Q N s (λ i ) 2 λ 3s i v 2 i + +∞ i=I min (N ) λ -s i -Q N s (λ i ) 2 λ 3s i v 2 i = Σ 1 (N ) + Σ 2 (N ).
The second sum Σ 2 (N ) can be bounded using eq. (4.86)

Σ 2 (N ) ε 2 (N ) 2 +∞ i=I min (N ) λ s i v 2 i . (4.89)
with

ε 2 (N ) = 2 sin(πs) π   π √ 2 e -(π 2 /2) √ N sin(πs) 1 -e -(π 2 /2) √ N + λ s-1 0 e 2(s-1) √ N 2(1 -s) (4.90) + e 2(1-s)/ √ N 2(1 -s) + 1 2s . So, lim N →+∞ ε 2 (N ) = sin(πs) πs(1 -s) . ( 4
λ s i v 2 i = 0. (4.93)
The first sum Σ 1 requires more efforts. By eq. (4.87) and using the convexity of the square function we have,

Σ 1 (N ) = I min (N )-1 i=1 λ -s i -Q N s (λ i ) 2 λ 3s i v 2 i 2 sin(πs) π 2 I min (N )-1 i=1   π √ 2 e -(π 2 /2) √ N sin(πs) 1 -e -(π 2 /2) √ N + λ σ-1 0 2(1 -s) e 2(s-1) √ N + 1 2s e -2sk(λ i )/ √ N 2 λ s i v 2 i 2 2 sin(πs) π 2 I min (N )-1 i=1     π √ 2 e -(π 2 /2) √ N sin(πs) 1 -e -(π 2 /2) √ N + λ σ-1 0 2(1 -s) e 2(s-1) √ N   2 + 1 2s e -2sk(λ i )/ √ N 2 λ s i v 2 i 8 sin(πs) 2 π 2   π √ 2 e -(π 2 /2) √ N sin(πs) 1 -e -(π 2 /2) √ N + λ σ-1 0 2(1 -s) e 2(s-1) √ N   2 I min (N )-1 i=1 λ s i v 2 i + 2 sin(πs) 2 π 2 s 2 I min -1 i=1 e -4sk(λ i )/ √ N λ s i v 2 i . (4.94) So, Σ 1 (N ) ε 1 (N ) 2 I min (N )-1 i=1 λ s i v 2 i + 2 sin(πs) 2 π 2 s 2 I min (N )-1 i=1 e -4sk(λ i )/ √ N λ s i v 2 i , (4.95)
with

ε 1 (N ) = 2 √ 2 sin(πs) π   π √ 2 e -(π 2 /2) √ N sin(πs) 1 -e -(π 2 /2) √ N + λ s-1 0 2(1 -s) e 2(s-1) √ N   , (4.96)
For the first term on the right-hand-side of eq. (4.95), by definition of ε 1 (N ) and I min (N ) and since u belongs to H s (Ω), we have

lim N →+∞ ε 1 (N ) = 0, (4.97)
and lim N →+∞ To do so, we introduce a new particular index defined by

I min (N )-1 i=1 λ s i v 2 i = v 2 H s . (4.98) So, lim N →+∞ ε 1 (N ) I min (N )-1 i=1 λ s i v 2 i = 0. ( 4
J min (N ) := min i ∈ 1, +∞ , λ i exp √ N . (4.101)
Since the eigenvalues of the Laplacian are increasing, we have J min (N ) I min (N ). Now, we can again split the sum from eq. (4.100) and get

2 sin(πs) 2 π 2 s 2 I min (N )-1 i=1 e -4sk(λ i )/ √ N λ s i v 2 i = 2 sin(πs) 2 π 2 s 2 J min (N )-1 i=1 e -4sk(λ i )/ √ N λ s i v 2 i (4.102) + 2 sin(πs) 2 π 2 s 2 I min (N )-1 i=J min (N ) e -4sk(λ i )/ √ N λ s i v 2 i .
We can show that the last sum tends to zero as N tends to +∞ follow,

λ s i v 2 i = 0. (4.105)
Consequently,

lim N →+∞ 2 sin(πs) 2 π 2 s 2 I min (N )-1 i=J min (N ) e -4sk(λ i )/ √ N λ s i v 2 i = 0. (4.106)
Now we bound the first sum in the right hand side of eq. (4.102). By definition of J min (N ) we have for each i in 1, J min (N ) -1 ,

1 2 ln(λ i ) 2 1 2 2 N 1 2 2 N, (4.107)
so since k(λ i ) = N -1 2 ln(λ i ) 2 , e -4sk(λ i )/ √ N e -4s √ N 1-1 2 2 = e -3s √ N -→ N →+∞ 0. (4.108)
Finally using eq. ( 4.108) in the first sum of the right-hand side of eq. (4.102),

2 sin(πs) 2 π 2 s 2 J min (N )-1 i=1 e -4sk(λ i )/ √ N λ s i v 2 i 2 sin(πs) 2 π 2 s 2 e -3s √ N J min (N )-1 i=1 λ s i v 2 i , (4.109)
and by eq. ( 4.104) and since v belongs to H s (Ω),

lim N →+∞ 2 sin(πs) 2 π 2 s 2 e -3s √ N J min (N )-1 i=1 λ s i v 2 i = lim N →+∞ 2 sin(πs) 2 π 2 s 2 e -3s √ N v 2 H s = 0. (4.110)
Consequently,

lim N →+∞ 2 sin(πs) 2 π 2 s 2 J min (N )-1 i=1 e -4sk(λ i )/ √ N λ s i v 2 i = 0. (4.111)
Combining eqs. (4.106) and (4.111) we obtain

lim N →+∞ Σ 1 (N ) = 0, (4.112)
thus,

lim N →+∞ v -v Q N s 2 H s lim N →+∞ +∞ i=1 λ -s i -Q N s (λ i ) 2 λ 3s i v 2 i lim N →+∞ Σ 1 (N ) + lim N →+∞ Σ 2 (N ) = 0.
This conclude the proof of Theorem 2 for the case σ = s.

The above result guarantees an exponential convergence rate of the sequence (v Q N s ) N to v in H s (Ω) only if v is more regular than H s (Ω) functions. When v belongs to H s (Ω) the convergence might be slower. In this case the bound we derive in the proof of Theorem 2 is given by

v -v Q N s 2 H s ε 1 (N ) 2 + 2 sin(πs) 2 π 2 s 2 e -3s √ N v 2 H s + 2 sin(πs) 2 π 2 s 2 + ε 2 (N ) 2 +∞ i=J min (N ) λ s i v 2 i , (4.113)
where ε 1 (N ) and ε 2 (N ) are respectively defined in eq. (4.96) and eq. (4.90) and where J min (N ) is defined in eq. (4.101). A similar result to Theorem 2 is given in [START_REF] Bonito | Numerical Approximation of Space-Time Fractional Parabolic Equations[END_REF]Theorem 7.1] when v = u solution to eq. (4.9) but it requires f ∈ L 2 (Ω) while Theorem 2 only requires f ∈ H s-2σ (Ω) for some σ ∈ [0, s) (especially, f belongs to negative order spectral fractional Sobolev spaces if σ > s/2) to ensure the exponential convergence rate.

Conjecture

The degenerated case σ = s in Lemma 1 limits the result in Theorem 2. In fact, the result of this theorem could be improved and the proof greatly simplified if the rational approximation Q N s was replaced by another rational approximation R N s satisfying a certain stronger convergence property (compared to eq. (4.15)).

We propose the following conjecture on the existence of such a rational approximation.

Conjecture 1. For s ∈ (0, 1) and λ 0 > 0, there exists a rational approximation scheme R N s such that:

R N s (λ) := C s N l=-N a l (1 + b l λ) -1 , ∀λ λ 0 , (4.114)
where the coefficients C s , {a l } N l=-N and {b l } N l=-N are positive real numbers and are independent of λ. Moreover, we assume

λ -s -R N s (λ) λ -s ε s (N ), ∀λ λ 0 , ∀N ∈ N * , (4.115)
with ε s independent of λ and lim N →+∞ ε s (N ) = 0 at an exponential rate.

If we assume Conjecture 1 and replace Q N s by R N s in eq. (4.73), we would obtain

v -v Q N s 2 H s +∞ i=1 ε s (N ) 2 λ s i v 2 i = ε s (N ) 2 v 2 H s . (4.116)
The convergence rate of the rational semi-discrete approximation v Q N s to the function v in the norm H s (Ω) would be preserved when v is exactly of H s (Ω) regularity, corresponding to a regularity of H -s for the data f when v = u.

Numerical results

We give numerical evidences showing that the rational approximation Q N s proposed in [START_REF] Bonito | Numerical approximation of fractional powers of elliptic operators[END_REF] does not satisfy eq. (4.115) in Conjecture 1. We compute the discrepancy between the functions λ → λ -s and λ → Q N s (λ) for a wide range of values of λ and up to large values of N , leading to very small outputs which can be polluted by rounding errors even with double-precision arithmetic. We avoid this difficulty by using the multi-arithmetic precision Python library gmpy2 [START_REF] Gladman | gmpy2 Python library[END_REF]. In fig. 4.1 we can see the plots of the function λ → λ s λ -s -Q N s (λ) . According to Conjecture 1 this function should be bounded uniformly by the quantity ε s (N ), converging to zero as N tends to +∞ and independent from λ. However, the plots in fig. 4.1 suggest that this is not the case for the rational approximation Q N s . Thus, Q N s does not seem to satisfy Conjecture 1.

Perspective

A future work will consist in investigating the truthfulness of Conjecture 1. As far as we know, no such rational approximation method satisfying eq. (4.115) is available in the literature. Although the property we are seeking is not plain uniform convergence, a direction that might worth considering is uniform rational approximation methods since many numerical evidences show their superiority over Q N s in terms of convergence [START_REF] Harizanov | Optimal solvers for linear systems with fractional powers of sparse SPD matrices[END_REF][START_REF] Harizanov | The Best Uniform Rational Approximation: Applications to Solving Equations Involving Fractional powers of Elliptic Operators[END_REF][START_REF] Hofreither | A unified view of some numerical methods for fractional diffusion[END_REF]. Another interesting extension to this work would be a posteriori error estimation. If we omit data approximation, eq. (4.113) in Theorem 2 is in fact a computable upper bound of the rational approximation error measured in the natural H s (Ω) norm when v = u (using eq. (4.6)). However, the function u Q N s is only a semi-discrete approximation to the solution u and a fully computable approximation to u can be obtained with the help of an additional approximation method, such as a finite element method. Then, in order to get a computable bound of the fully discrete approximation, a posteriori estimation of the finite element discretization error is needed. To balance the discretization errors coming from the rational approximation method and from the finite element method would allow to derive efficient numerical methods to solve fractional Laplacian equations.

To our knowledge only very few works on a posteriori error estimation applied to fractional problems are available in the literature [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF][START_REF] Antil | Fractional operators with inhomogeneous boundary conditions: analysis, control, and discretization[END_REF][START_REF] Bonito | Numerical methods for fractional diffusion[END_REF][START_REF] Nochetto | A PDE Approach to Fractional Diffusion in General Domains: A Priori Error Analysis[END_REF]. Except in [START_REF] Ainsworth | Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver[END_REF], a posteriori error estimators are all based on the Caffarelli-Silvestre extension allowing to reformulate the fractional problem eq. (4.4) into a non-fractional problem on Ω × R, which complicates its implementation when Ω is a three-dimensional domain.

Introduction

Human meniscus plays a key role in the functioning of the knee joint fig. 5.1 a). This tissue has a number of functions such as: load bearing (about 45 -75% of the total load on the joint), joint stability and lubrication [START_REF] Kurosawa | Load-bearing mode of the knee joint: physical behavior of the knee joint with or without menisci[END_REF][START_REF] Shrive | Load-bearing in the knee joint[END_REF]. Degenerative processes of the meniscus either from injuries or aging affect approximately 35% of the population [START_REF] Sihvonen | and for the Finnish Degenerative Meniscal Lesion Study Group. Mechanical Symptoms and Arthroscopic Partial Meniscectomy in Patients With Degenerative Meniscus Tear: A Secondary Analysis of a Randomized TrialMechanical Symptoms and Arthroscopic Partial Meniscectomy[END_REF]. When required, the most common surgical procedure is currently total/partial meniscectomy (i.e removal of the damaged tissue). Therefore, partial/total meniscal replacements are expected to help avoid articular cartilage degeneration. Currently the clinical outcomes of these implants are not ideal [START_REF] Van Kampen | The knee joint in sports medicine[END_REF] due to the fact that they do not mimic the structureproperty relationships of the tissue as these are still not well understood [START_REF] Fox | The basic science of human knee menisci: structure, composition, and function[END_REF]. Among many biologic tissues, the meniscal tissue is composed of porous solid matrices -mainly collagen-with fluid filling the pores [START_REF] Agustoni | High resolution microcomputed tomography reveals a network of collagen channels in the body region of 113 the knee meniscus[END_REF][START_REF] Bonomo | A procedure for slicing and characterizing soft heterogeneous and irregular-shaped tissue[END_REF][START_REF] Maritz | The functionally grading elastic and viscoelastic properties of the body region of the knee meniscus[END_REF][START_REF] Vetri | Advanced microscopy analysis of the micro-nanoscale architecture of human menisci[END_REF]. The overall mechanical behavior of this type of tissues depends not only on the solid matrix deformation, but also on the movement of the fluid in and out of the collagen channels [START_REF] Agustoni | High resolution microcomputed tomography reveals a network of collagen channels in the body region of 113 the knee meniscus[END_REF][START_REF] Maritz | The functionally grading elastic and viscoelastic properties of the body region of the knee meniscus[END_REF] during the deformation. Meniscal tissue exhibit a non-uniform and anisotropic porosity, which is related to graded material properties. Such graded material properties are fundamental for the correct functioning of this tissue [START_REF] Maritz | The functionally grading elastic and viscoelastic properties of the body region of the knee meniscus[END_REF]. The investigation of the material properties, coupled with the quantification of architectural parameters such as porosity and channels interconnectivity, will enable the design of artificial cellular structure, which can resemble the native behavior of the tissue (see for instance [START_REF] Nazir | A state-of-the-art review on types, design, optimization, and additive manufacturing of cellular structures[END_REF]). More in depth, advanced microscopy investigations highlighted that collagen bundles form the wall of channels which can be observed both at the macroscale and at the microscale [START_REF] Agustoni | High resolution microcomputed tomography reveals a network of collagen channels in the body region of 113 the knee meniscus[END_REF][START_REF] Vetri | Advanced microscopy analysis of the micro-nanoscale architecture of human menisci[END_REF]. Fluid is able to flow inside these channels in response to physiological loading. Understanding and appropriately modeling the fluid flow behavior in the different portions of the meniscal tissue, considering a range of various loading conditions, is essential to gain insight on the biomechanical function of this tissue.

To date, limited information is available on region-specific and anisotropic permeability in the meniscus, which is essential to understand the fluid flow evolution and its relationship with the internal architecture.

Poromechanics experimental tests, such as confined compression tests -including both stress relaxation and creep-are currently used in order to characterize material parameters, such as elastic modulus at equilibrium and permeability. In order to identify these two parameters the main models used are based on the biphasic and poroelastic theories.

Even though these theories developed from different roots, it can be shown that they are basically equivalent [START_REF] Bowen | Compressible porous media models by use of the theory of mixtures[END_REF].

In this study, we focused our attention on the evaluation of the permeability in the meniscus and in particular on its variation in the different regions of the tissue and directions (See fig. 5.1 b). Moreover, we showed that the biphasic model does not provide a good fit with the experimental curves, therefore we proposed a poroelastic model in which the pore-pressure diffusion equation is derived by adopting a modified version of Darcy's law involving fractional derivatives [START_REF] Barrera | A unified modelling and simulation for coupled anomalous transport in porous media and its finite element implementation[END_REF]. It has been shown that in a high porosity medium there is a departure from the Darcy's law as the inertia (velocity-squared term), thermal dispersion, convective (development term) and boundary (no-slip condition) effects not included in the Darcy's model may play a significant role [START_REF] Jang | Variable porosity and thermal dispersion effects on vortex instability of a horizontal natural convection flow in a saturated porous medium[END_REF]. Moreover, in the last few decades experimental evidence of anomalous diffusive phenomena, i.e. not following a Darcian behavior, has grown [START_REF] Abd | Neutron radiography study of water absorption in porous building materials: anomalous diffusion analysis[END_REF][START_REF] De Azevedo | Water ingress in y-type zeolite: Anomalous moisture-dependent transport diffusivity[END_REF][START_REF] Iaffaldano | Experimental and theoretical memory diffusion of water in sand[END_REF][START_REF] Küntz | Experimental evidence and theoretical analysis of anomalous diffusion during water infiltration in porous building materials[END_REF][START_REF] Ramos | Anomalous diffusion during water absorption in porous building materials -experimental evidence[END_REF]. This is mainly due to the fact that the permeability, hence the rate of fluid flow, is not a constant quantity. Variations in permeability occur, for example, when the fluid flow impacts the geometry or the micro-structural features, such as the configuration of the pores. For example, experiments on water flow in building materials highlighted that the permeability changes during the flowing process as a result of the microstructural rearrangement of grains/pores. Iaffaldano at al. [START_REF] Iaffaldano | Experimental and theoretical memory diffusion of water in sand[END_REF] hinted that during compaction of sand, permeability might decrease due to the fact that the fluid carries solid particles which then close some of the pores. Essentially, the configuration of the medium, in particular the ratio between closed/open pores, changes during the process. Fluid might be trapped in the medium leading to a slower fluid flow rate. On the contrary, if during the fluid diffusion process some of the pores open creating conductive microchannels, permeability might increase. Therefore, fluid can be transported for large distance in a reduced time determining a faster diffusion process. Modeling the anomalous fluid diffusion process is one of the key points when dealing with poromechanics of biological tissues. Therefore, it is fundamental to develop a theory, which can incorporate the change in microstructural features (for instance, the interaction between fluid particle and open pores/channels) during the transport process. Recently the stochastic Continuous Time Random Walk (CTRW) framework was proposed for this purpose [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF]. However a deterministic derivation, dual to the CTRW, can be represented by the introduction of a modified version of Darcy's law involving linear fractional operators [START_REF] Caputo | Diffusion in porous layers with memory[END_REF][START_REF] Płociniczak | Analytical studies of a time-fractional porous medium equation. derivation, approximation and applications[END_REF]. The goal of this chapter is to present the results of experimental confined compression tests performed on samples extracted from the three portions of the human knee meniscus (posterior, central and posterior). During confined compression, fluid flows through the collagen channels with a rate depending on the permeability of the tissue itself (fig. 5.1 c). We aim at: (1) extracting the diagonal terms of the permeability tensor; (2) studying, analyzing and modeling the observed decay of fluid flow during the test, which is not captured by the classical Darcy's law. Hence, we propose a generalization of such law which involves time derivatives of non-integer order appropriate to model the fluid flow in the meniscal tissue; (3) proposing a small deformation fractional poroelastic model for the human meniscal tissue during confined compression. The model enables us to identify the two parameters involved in the fractional pore pressure diffusion equation, namely the permeability and the order of the (fractional) derivative. In this study, we do not consider coupling between the flow in different directions, i.e. assume that the permeability tensor is diagonal. The structure of the chapter is as follows: we introduce the rationale behind a fractional Darcy's law, we then summarize the main equations of both biphasic and fractional poroelastic theories; we then present the confined compression poromechanics experimental tests and discuss the material parameters (fractional permeability and order of the fractional derivative) we recover through the fittings.

Fluid flow in complex porous media: time-fractional Darcy's law

Fluid flow in porous media is commonly modeled by Darcy's law. This relation was observed experimentally by Darcy [START_REF] Darcy | Les Fontaines publiques de la ville de Dijon[END_REF] and derived in 1986 from the Navier Stokes equations using homogenization theory [START_REF] Whitaker | Flow in porous media i: A theoretical derivation of darcy's law[END_REF]. Darcy's law states that the instantaneous flow rate through a homogeneously permeable porous medium of permeability k is proportional to the dynamic viscosity of the fluid and the pressure drop over a given distance. The total discharge, J f (units L 3 /T , where L and T indicate units of length and time respectively) equates the product of the intrinsic permeability, k (L 2 ) with the cross-sectional area flow, A (L 2 ) and the total pressure drop p out -p in (F /L 2 ), all divided by the dynamic viscosity, µ (F T /L 2 ) and the length over which the pressure drop is taking place :

J f = - k • A • (p out -p in ) µ
More generally, Darcy's law states that the discharge per unit area, j f = J f /A is proportional to the pressure gradient, the intrinsic permeability and inversely proportional to the dynamic viscosity:

j f = - k µ ∇p.
The derivation of Darcy's law from Navier Stokes equations assumes a creeping, laminar, stationary and incompressible flow of density ρ and velocity (u i ) 1≤i≤3 . Incompressibility implies: ( Dρu i Dt = 0), which leads to Stokes' equation, in the presence of gravity g i :

µ∇ 2 u i + ρg i -p ,i = 0.
Assuming that the viscous resisting force varies linearly with the velocity, introducing the porosity ϕ, and the second order permeability tensor k = k ij , a simple derivation leads to:

-(k ij ) -1 µϕu j + ρg i -p ,i = 0.
The discharge per unit area in direction n can then be written as follows:

j f n = - k ni µ (p ,i -ρg i ) . j f = - k µ (∇p -ρg) .
If the resistance to fluid flow offered by the pores varies in space, the components of the permeability tensor k ij vary in space. k ij is a symmetric (Onsager reciprocal relations), positive definite (because the flow component parallel to the pressure drop occurs in the same direction as the pressure drop) matrix. k ij may be isotropic, in which case it is diagonal and all diagonal entries are identical: k ij = kδ ij . In general, the permeability tensor is anisotropic, and may also not be diagonal. In all cases, the permeability tensor can be diagonalized, as it is symmetric positive definite. As previously underlined, when dealing with soft tissue and complex porous media, permeability changes during the deformation process, therefore, it is important to develop a theory, such as the CTRW, which incorporates the change in microstructural features during the transport process. A modified version of Darcy's law involving linear fractional operators can be seen as the deterministic version of CTRW [START_REF] Caputo | Diffusion in porous layers with memory[END_REF][START_REF] Płociniczak | Analytical studies of a time-fractional porous medium equation. derivation, approximation and applications[END_REF]:

j f = - k µ D β 0 (∇p) (5.1)
In the following k µ will be indicated as λ β . D β a indicates the Caputo's fractional derivative of order β and origin a of ∇p which is defined below:

(D β a ∇p)(t) = 1 Γ(n -β) t a (t -τ ) n-β-1 ∂ n ∂τ n ∇p(τ )dτ (5.2)
The expression is valid for n -1 < β < n, and Γ is the Euler's Gamma function. In the cases considered in this chapter the origin a = 0. The fractional derivative method offers the possibility to model, with reduced number of parameters, all of the anomalous diffusion behaviors by changing the order of the derivative. The main drawback is that it is difficult to link the order of the derivatives with the microstructural features. Within this frame, in this chapter we identify a simple mathematical model able to describe fluid flow in the human meniscus.

Biphasic and Linear fractional Poroelastic models

The assumptions at the basis of both the poroelasticity and biphasic theories are the following:

• the solid phase is incompressible, linear elastic, subject to infinitesimal strain, homogeneous and isotropic (i.e., the material parameters of the solid phase do not depend on the orientation nor the position in the sample), non-dissipative;

• the fluid phase is incompressible, it flows slowly through the pores, it is homogeneous, isotropic and non-dissipative and there is no fluid source;

• the only dissipation comes from the frictional drag due to the relative velocities of the two phases;

• the absence of external body forces (other than those explicitly mentioned for the confined compression tests);

• the isotropy and homogeneity of the permeability tensor, which is then reduced to the scalar k representing the averaged intrinsic permeability of the sample;

• the permeability k is a constant parameter (time-independent) in the biphasic model and it is a time-dependent quantity in the fractional poroelastic model (due to the time-fractional Darcy's law described in eq. (5.1)).

Biphasic model -consolidation problem

For the purposes of this study, we restrict our formulation to small strain theory and consider the linear biphasic model from [START_REF] Mow | Biphasic Creep and Stress Relaxation of Articular Cartilage in Compression: Theory and Experiments[END_REF]. We assume the solid matrix is incompressible, linear elastic, isotropic, homogeneous and non-dissipative, whereas the interstitial fluid is incompressible and non-dissipative. The only dissipation comes from the frictional drag due to the relative velocities of the two phases. We denote λ 0 = k/µ the averaged axial permeability of the sample and H A the aggregate modulus. Under these assumptions and the additional assumptions coming from the confined compression tests setting, the biphasic theory leads to the following unidimensional boundary value problem for the vertical displacement u s z of the solid phase

∂ 2 u s z ∂z 2 = 1 H A λ 0 ∂u s z ∂t in (0, h) × (0, T ), (5.3)
the boundary condition at z = 0 where the solid skeleton is fixed

u s z (z = 0, t) = 0 on (0, T ), (5.4)
the initial condition u s z (z, t = 0) = 0 on (0, h),

and finally the boundary condition at z = h which will be different for the creep test and for the stress relaxation test [START_REF] Soltz | Experimental verification and theoretical prediction of cartilage interstitial fluid pressurization at an impermeable contact interface in confined compression[END_REF]. For the creep test the condition is

∂u s z ∂z |z=h = - P A H A on (0, T ), (5.6)
where P A is the applied compressive stress and for the stress relaxation test the condition is

u(z = h, t) = -V 0 t if 0 t < t 0 , -V 0 t 0 if t 0 t T, (5.7)
where V 0 and t 0 are input data from the stress relaxation test.

If we denote c n = (-1)

n-1 2
2h nπ 2 the solution to eqs. (5.3) to (5.6) in the case of the creep test reads

u s z (z, t) = P A H A   -z + 2 h ∞ n=1,3 c n exp - H A λ 0 n 2 π 2 t 4h 2 sin nπz 2h   , (5.8)
In the stress relaxation test case, if we denote ρ n = nπ h 2 , the solution to eqs. (5.3) to (5.5) and (5.7) [START_REF] Soltz | Experimental verification and theoretical prediction of cartilage interstitial fluid pressurization at an impermeable contact interface in confined compression[END_REF] is given by

u s z (z, t) =          - V0tz h - 2V0 HAλ0h ∞ n=1 (-1) n ρ 3/2 n 1 -e -H A ρnλ 0 t sin ( √ ρnz) if 0 t < t0, - V0t0 h - 2V0 HAλ0h ∞ n=1 (-1) n ρ 3/2 n e -H A ρnλ 0 t e -H A ρnλ 0 t -1 sin ( √ ρnz) if t0 t T.
(5.9)

The stress relaxation response σ t is given by H A ∂u ∂z |z=h , so:

σ t =            - V 0 H A t h + 2V 0 λ 0 h ∞ n=1 ρ -1 n 1 -e -H A ρnλ 0 t if 0 t < t 0 , 2V 0 λ 0 h ∞ n=1 ρ -1 n e -H A ρnλ 0 t e -H A ρnλ 0 t -1 if t 0 t T.
(5.10)

Linear fractional poroelastic model: fractional consolidation problem

The where λ = K -2 3G is the Lamé constant and G, K are the shear and bulk modulus respectively, p is the pore pressure, α is the Biot coefficient. In order to solve this one, we need an additional equation, which is given by the pore pressure diffusion equation derived considering a time fractional Darcy's law in eq. (5.1) [START_REF] Barrera | A unified modelling and simulation for coupled anomalous transport in porous media and its finite element implementation[END_REF] :

∂p ∂t = KB α λ β D β 0 ∇ 2 p -B ∂σ H ∂t (5.12)
Where D β 0 indicates the Caputo's time fractional derivative [START_REF] Podlubny | Fractional differential equations : an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF]. It is important to note that the classical (non-fractional) pore pressure diffusion equation is recovered in case of β = 0. In this case λ 0 = k µ with dimension of

[L] 4 [F ][T ] .
In the case of β = 0 note that λ β has dimension of [L] 4

[F ][T ] 1-β . The pore pressure diffusion equation can also be written in terms of strains [START_REF] Barrera | A unified modelling and simulation for coupled anomalous transport in porous media and its finite element implementation[END_REF]:

∂p ∂t = K u -K α 2 λ β D β 0 ∇ 2 p - K u -K α ∂ε d ∂t (5.13)
The consolidation problem is modeled through the 1D uniaxial strain poroelastic problem [START_REF] Barrera | A unified modelling and simulation for coupled anomalous transport in porous media and its finite element implementation[END_REF]. Equation (5.11) adapted for the 1D case in which the only non-zero component of strain is ε zz , we obtain:

K + 4G 3
∂ε zz ∂z -α ∂p ∂z = 0 (5.14)
Equation (5.14) is then coupled with the pore pressure diffusion equation in eq. (5.13) in order to obtain the following pore pressure diffusion equation:

∂p ∂t = λD β 0 ∂ 2 p ∂z 2 (5.15)
where λ = λ β (4G+3K)(Ku-K) α 2 (4G+3Ku) . We indicate with the symbol λ 0 = (4G+3K)(Ku-K) α 2 (4G+3Ku) . The boundary value problem is given [START_REF] Barrera | A unified modelling and simulation for coupled anomalous transport in porous media and its finite element implementation[END_REF]:

∂p ∂z z=0 = 0, (5.16)
p(z = h, t) = 0, (5.17)

u(z = h, t) = 0.

(5.18)

A constant compressive stress in the z direction is applied to the cylinder at z = 0:

σ zz (0, t) = -P A (5.19)
where -P A is the applied compressive stress. The initial pore pressure is derived for undrained conditions i.e:

p(z, 0) = P A 3(K u -K) α (4G + 3K u ) (5.20)
The analytical solution in terms of pore pressure reads as follows:

p(z, t) = P A γ ∞ n=1,3 E 1-β,1 - n 2 π 2 λt 1-β 4h 2 c n cos nπz 2h (5.21)
where:

γ = 3 (K u -K) α (4G + 3K u ) (5.22) c n = (-1) n-1 2 2h nπ 2 (5.23)
where E 1-β,1 is the Mittag-Leffler function. In the case of β = 0 the solution is identical to the classical Terzaghi's solution in which E 1-β,1 = exp.

In the case of β = 0, it is possible to derive the following displacement analytical solution [START_REF] Detournay | Fundamentals of Poroelasticity[END_REF]:

u z (z, t) = P A λ 0   -z + 2 h ∞ n=1,3 c n exp - λn 2 π 2 t 4h 2 sin nπz 2h   , (5.24)
where c n = (-1)

n-1 2 2h nπ 2 .

Correspondence of parameters between the biphasic and the linear fractional poroelastic model

The

((4G+3K)(Ku-K)) α 2 (4G+3Ku) m 3 • s • kg -1 -incompressible case 1 3 (4G + 3K) λ β fractional averaged axial permeability m 3 • s 1+β • kg -1 β
fractional power Unitless depends instead on three parameters: the diffusion coefficient λ β , the fractional power β and the ratio k/µ. These parameters depend on the material parameters of the constituent phases or have to be determined, for example using numerical experiments curves fitting. The biphasic and linear fractional Biot models are equivalent under the above assumptions when β = 0 (i.e. in the case of a non-fractional Biot's model) and when H A = λ 0 . Equations (5.8) and (5.24) show that biphasic and Biot models are equivalent when H A = λ 0 . Similarly, the displacement solution from the Biot's model in the case of the stress relaxation test can be obtained by replacing H A by λ 0 in eq. (5.9). Tables 5.1 and 5.2 summarize the relationships between the parameters for each models.

Materials and methods

Poromechanics tests -Confined compression

Menisci were harvested from patients (age 65-76, mean 72, standard deviation 4) undergoing total knee arthroplasty (ethical approval EM 249-2018 21/2017/Sper/IOR EM2, Rizzoli Orthopaedic Institute, Bologna, Italy). Samples labeled as "degraded" by gross investigation of the surgeon were discarded. Three lateral menisci and three medial menisci were collected and stored at -20 C [START_REF] Forkel | Biomechanical and viscoelastic properties of different posterior meniscal root fixation techniques[END_REF][START_REF] Chia | Compressive moduli of the human medial meniscus in the axial and radial directions at equilibrium and at a physiological strain rate[END_REF]. The day of the test, each meniscus was thawed in a phosphate-buffered saline (PBS) bath at room temperature for about thirty minutes [START_REF] Chia | Compressive moduli of the human medial meniscus in the axial and radial directions at equilibrium and at a physiological strain rate[END_REF]. Then, cylindrical samples (diameter of 3mm, height 3-4 mm) were extracted from the central body, anterior and posterior region along three reference directions, i.e. vertical, radial and circumferential (fig. 5.1 b), following a dedicated procedure [START_REF]Mach-1 -confined compression of a cartilage disk[END_REF], adapted to the meniscal configuration. A total of 18 cylindrical samples were grouped considering the the harvesting region (body, anterior and posterior) and the direction (vertical, radial, circumferential) and then tested. The testing protocol was implemented on a multi-axis mechanical tester (Mach-1, Biomomentum Inc., Canada) in a confined compression configuration fig. 5.1c); with this specific setup, we were confident that the fluid could flow only from the base of the cylinder. Insertion in the confining chamber, thickness measurement and removal of the meniscal sample followed a dedicated procedure [START_REF]Mach-1 -confined compression of a cartilage disk[END_REF]. Concerning confined compression, the setup involved a confining chamber with an inner diameter equal to that of the tool used during the extraction of cylindrical samples, i.e. 3 mm. The bottom of the chamber consisted of a porous-permeable platen, while the top allowed the insertion of the piston. Both these components (Biomomentum Inc., Canada) were manufactured to allow the leakage of the fluid only through the porous platen. Moreover, an additional control was provided by visual inspection of the confining chamber transparent wall, which allowed to check any fluid flow towards the piston (fig. 5.2 a). We measured the amount of fluid flowing out each sample as result of the applied compression loading by monitoring the weight of the samples throughout the test.

Testing protocols

Relaxation tests. For the relaxation tests, we followed the procedure for permeability and aggregate modulus H A analysis, as recommended by the manufacturer of the testing machine [START_REF]Mach-1 -confined compression of a cartilage disk[END_REF], [START_REF]Mach-1 -extraction of mechanical parameters following confined compression[END_REF]. More in detail, the test sequence considers a pre-compression with ramp amplitude 10% of the sample thickness (h) and ramp velocity 0.3%h/s, followed by five stress-relaxations with incremental ramp amplitude 2%h and ramp velocity 0.3%h/ s.

Creep tests

The testing protocol -specifically the loading amplitude and the detection of the fluid flow out of the sample -was designed and optimized for this specific study.

In particular, three testing phases were implemented for each meniscal sample. The first phase, developed according to [START_REF]Mach-1 -confined compression of a cartilage disk[END_REF], consisted in five separated repetitions of confined compression, 75 s each, thus to realize 450 s of total creep. More in detail, before and after each compression, sample was removed from the confining chamber, weighed by a microbalance (Tecnopound, Ravenna, Italy) and repositioned inside the chamber with the Bulk modulus K=1.6x10 5 Pa shear modulus G=76923 Pa (E = 0.2 • 10 6 Pa, ν = 0.3) Skempton coefficient B=0.88 Biot coefficient α=0.65 Undrained bulk modulus K u =K/1α B Table 5.3: Material parameters previous vertical orientation. The second phase involved the resting of the sample in PBS until its height returned to the pre-loading value, thus to recover the loading history. In the final phase, i.e. third one, a single step of confined compression creep was applied to the sample for 450 s. Similarly to what we reported for the phase one, also in this case the weight of the sample was measured before and after the compression. It is crucial to emphasize that this third phase served to calibrate the weights measured during the phase one, when the multiple removals of the sample from the confining chamber and the separated compressive ramps could have affected the creep. Stage velocity was 0.3 % h, while the load target was 0.5 N, corresponding to a stress of about 0.07 MPa, which is within the range of physiological values for human menisci [START_REF] Seitz | Stress-relaxation response of human menisci under confined compression conditions[END_REF]. The decrease in weight of the sample has been monitored throughout the test and it has been correlated to the amount of fluid discharged by the samples, hence it served to reveal information regarding the rate of the fluid flow. In particular, the weight of the sample over time W (t) is related to the initial weight of the sample W 0 and the fluid flux:

W (t) = W 0 - t 0 j f • w s Adt (5.25)
With j f given by eq. (5.1), in which the gradient of the pressure is given by differentiating the expression of the pressure in eq. (5.21), w s being the specific weight and A the crosssectional area. Substituting eq. (5.1) and eq. (5.21) into eq. ( 5.25) we obtain the following relation:

W (t) = W 0 -P A γ 2 h ∞ n=1,3 λ β t 1-β E 1-β,2-β - π 2 λt 1-β 4h 2 (5.26)
In eq. (5.26) n=1 has been retained since no improvements are observed by adding more terms to the series. This means that one term is enough to describe the physics of this kind of test, but in more realistic conditions it is expected that more terms of the summation are needed to accurately describe the solution. 2 contains the parameters H A and permeability k of all the 5-steps.

Poromechanics test results and fittings

Relaxation tests results -biphasic model

By implementing the standard procedure (i.e., confined compression, stress relaxation, no weight measurements and classic biphasic theory) -which was designed for cartilage, but it is here adopted for meniscal cylinders -the revealed mechanical parameters presented a large variability (fig. 5.3). Moreover, the expected trend of decreasing permeability with increasing strain was visible only for one trial ("TK11-CB-Vert" in table 5.4). This fact, coupled with high values of fitting root mean squared error (RMSE), underlines the difficulties of this approach in characterizing meniscal tissue.

Creep tests results -fractional poroelastic model

Focusing on confined compression tests, we specifically measured the decrease in weight of the cylinder during the tests and correlate this data to the amount of fluid discharged by the samples j f during the test. We now illustrate the best-fitting procedure performed

Ramp amplitude Sample

Step (% height) between the proposed model in eq. (5.26) and the experimental results as shown in fig. 5.4 a-f for a few samples taken from one of the three lateral menisci. In particular, two free parameters are considered: β and λ β . We perform the best-fitting procedure for a time window extended from t = 0 to 450s and γ = 0.695 has been fixed by considering specific mechanical parameters which are summarized in table 5.3. The values of the parameters obtained through the fittings are reported in table 5.5.

H A (MPa) k (m 4 /s × N) RMSE 1

Discussion

In this study, we assumed that the fluid flow is ruled by a modified version of Darcy's law (eq. (5.1)). According to eq. (5.1), the fluid flow is not steady as modeled by the classical Darcy's law. Instead the fluid flow rate evolves with time, more specifically it evolves with a fractional time derivative (of order β). Equation (5.1) also implies that the permeability is anomalous in the sense that its units are function of the order of the derivative λ β = [ L 4 F T 1-β ]. The classical Darcy's law is recovered for β = 0 in eq. (5.1). In order to estimate how closely both models (classical and fractional Darcy's law) fit the experimen- tal data we compared the RMSE values. We observe in 12/18 fitted experimental tests reported in table 5.5 the RMSE given by the fractional model is lower by 10% with respect to the classical model. Although this might seem not enough to justify a new theory, it is important to note that this is the first attempt to use this type of model for interpreting the rate of fluid flow inside a biological tissue. We are confident that a more precise experimental set up will give more emphasis of the benefit of having a generalized fractional model which allow to incorporate the classical one by simply setting β=0. Furthermore, it would be interesting to study the evolution of fluid flow inside the meniscus and how it varies both spatially (within the posterior/central/anterior portions) of the tissue and directionally (in the vertical/radial/circumferential directions) considering the parameters in table 5.5. A first analysis of the parameters highlights that the values of the order of the fractional derivative β and the anomalous permeability It can be noted that, although the values of β in the three regions might not significantly diverge from zero, they indeed affect the evolution of the fluid flow rate. In order to investigate the role of the order of the fractional time derivative β, we specifically realized a computational simulation in which we applied a constant gradient of pressure ∇p and we calculated the fluid flow rate by applying the fractional Darcy's law relation as given in eq. (5.1). it is possible to observe that the higher the value of β the faster the decrease in the fluid flow rate is. Furthermore, we analyzed the evolution of the fluid flow rate in the central body, posterior and anterior horns using the mean values of the parameters β, λ β above. Figure 5.6b shows the response of the normalized flux in the three regions. It can be noted that as the value of the anomalous permeability λ β in the central body is about four times higher than in the anterior horn, the value of flux is higher in the central body and lowest in the anterior region. Moreover, as the value of β is larger compared to the anterior and posterior horns, the fluid flux is faster (i.e. decrease in fluid flow rate) in the central body. Ns (1-β) , respectively. Figure 5.6c pictures the evolution of the normalized fluid flow rate in the three directions. It can be noted that the values of β are very close for all the three directions, differing within a few percents. Therefore the decrease of flux in time (i.e., the flux velocity) can be considered as almost identical in the three main directions. The value of β rules the time evolution of the pore pressure diffusion in eq. (5.21). Increasing the value of β implies a faster pore pressure diffusion with time. At the beginning of the test, when the fluid saturates the pores, the pore pressure carries most of the load. As the test continues, the fluid flows out of the specimen. The pore pressure decreases and hence the solid structure starts deforming. The rate at which the fluid flows and the pore pressure decrease and the solid structure deformation increases is ruled by the value of the fractional derivative β. A higher value of β implies a faster pore pressure diffusion and hence a faster solid deformation. Figure 5.7 shows the evolution of the adimensional pore pressure (p = p(z, t)/P A in eq. (5.21)) throughout the length of the sample along the z axis for different values of β = 0, 0.1, 0.5 at time 15s. It can be seen that the pore pressure drops at a faster rate with increasing β.

Conclusions

This work focuses on the understanding of the evolution of the fluid flow inside the human meniscal tissue. We showed through µCT scans of the meniscus that the porosity significantly varies spatially within small portions of the tissue. This leads to a functionally graded permeability across the meniscal tissue. We wanted to address a question regarding how to appropriately model the fluid flow inside the meniscus. We specifically performed confined compression tests on samples extracted from three portions (posterior/central/anterior) of the structure and in three directions (vertical / radial / circumferential) in order to measure experimentally the anisotropic permeability of the human meniscal tissue. We then correlated the weight loss of the sample with the evolution of the fluid discharged upon compression of the tissue. Results show that the weight loss of the sample is well described by a three-parameters equation derived from a fully coupled poroelastic model in which the fluid flux evolution is ruled by a generalized Darcy's law involving fractional operators such as derivatives of non-integer order. We obtained the anisotropic parameters needed to characterize the fluid flow evolution inside the different parts of the meniscus, i.e. the "anomalous" permeability as well as the order of the derivative. We noticed that the parameters of this porous medium equation are functionally graded in space. Our results preliminary suggest that the fluid flow in the central part is faster than in the posterior and anterior regions. Moreover, we noted that the flux is higher in the circumferential direction of the central body compared to the radial and vertical directions. However, the decrease in flux over time (i.e. flux velocity) can be considered almost identical in the three directions. We believe that this work is a fist attempt and a pioneering study which experimentally investigates the anomalous behavior of the meniscal tissue. Furthermore, the presented approach can be easily adapted to study other types of biological tissues. On the other hand, there are several limitation of this study, as hereinafter summarized:

• additional experimental tests (unconfined compression, confined compression with control of the pressure) need to be performed in order to fully characterize the single poromechanics parameters appearing in λ β and λ;

• a wider sample size and additional information is required to perform a reliable statistical analysis (so far we tested one cylinder per portion of meniscus in the vertical, radial and circumferential directions).

• This work focuses on the 1D poromechanics behavior in the three directions, however in order to build the full anomalous permeability tensor, the coupled behavior should be assessed in depth. Nevertheless, this study requires a more complex experimental set up.

Perspective

No matter how accurate scientific models and theories become, the quantification of uncertainties will remain of crucial importance. In this dissertation we have contributed to improve and extend some techniques of error estimation in the context of finite element methods with a particular emphasis on fractional Laplacian equations. Despite the extensive literature on the Bank-Weiser estimator, some open questions remain. We can cite particularly the proof of Property i) without the use of the saturation assumption for the general definition of the Bank-Weiser estimator. Especially for finite elements of degree higher than one. The dependence of the Bank-Weiser estimator performance on the choice of the spaces used in the estimator definition would also be worth to investigate further. The derivation and study of Bank-Weiser-like estimators for finite elements other than Lagrange elements is also to be explored. Finally, a convergence proof of adaptive refinement algorithms steered by the Bank-Weiser estimator still remains to be derived.

In this manuscript we have introduced the first FEniCS library dedicated to a posteriori error estimation: FEniCS-Error-Estimation. The library is mostly focused on the implementation of the Bank-Weiser estimator, only the simplest versions of two other kinds of estimators (explicit residual and Zienkiewicz-Zhu) being implemented in the library. In the future, it would be interesting to add new a posteriori error estimation techniques to this library, especially the methods giving proper boundary conditions to the local problems in order to derive guaranteed implicit residual estimators. Flux equilibration estimators would be another interesting and useful addition to the library. Beside other a posteriori error estimation techniques, new test cases would be worth to implement such as time-dependent equations, non-linear equations or equations based on non-local and fractional operators. For the moment only Lagrange finite elements have been implemented. Thus, adding a posteriori error estimators for other kinds of elements would also be an interesting future direction.

In Chapter 3 we have introduced a novel a posteriori estimator for the L 2 error induced by a finite element discretization of spectral fractional Laplacian operators, based on a particular rational approximation method. In the future it would worth considering the generalization to other types of rational approximations (following e.g. [START_REF] Hofreither | A unified view of some numerical methods for fractional diffusion[END_REF]). To adapt the estimator to integral fractional Laplacian equations would be another interesting generalization. This could be achieved, for example, following [START_REF] Bonito | Numerical methods for fractional diffusion[END_REF]Equation 1.13]. The mathematical justification of our estimator is missing and would be worth to investigate. It is well-known that the justification of a posteriori estimators of the error is often more demanding in the L 2 -norm than in the energy norm (the spectral fractional Sobolev norm in our case) since it requires the use of the Aubin-Nitsche argument which demands a sufficient regularity on the solution (see e.g. [START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF]). So, another very interesting direction would be to adapt this method to the estimation of the spectral fractional norm of the error. Finally, from a computational perspective it would be interesting to control the error on the rational approximation scheme in order to balance the two sources of error: the rational approximation one and the finite element discretization one. In addition, our method allows to estimate the error committed on each parametric problem of the rational decomposition and therefore allows to derive a specific adaptive strategy for each one of these parametric problems. This would reduce the computational cost of the overall method e.g. by calibrating a mesh for each parametric problem.

An interesting extension to Chapter 4 would be to combine the error estimation of the rational approximation method with an a posteriori error estimation of a finite element method in order to derive an efficient adaptive mesh refinement strategy for the discretization of the spectral fractional Laplacian.

In Chapter 5 we proposed the first experimental investigation of the anomalous behavior of the meniscal tissue. Additional experimental tests need to be performed to fully determine the poromechanics parameters of the model. In particular, additional information is required to make a reliable statistical analysis. Our work only treats one-dimensional experiments (in the three directions). The characterization of the full permeability tensor requires a three-dimensional study which involves a more complex uncoupled poro-elastic model and necessitates a more complex experimental set up. We only tackled the nonlocality in time (memory) of the tissue, an interesting extension would be to investigate non-locality in space via e.g. fractional Laplacian equations.

The class of residual estimators, the explicit residual estimator is part of, have been introduced for the first time in [START_REF] Babuška | A-posteriori error estimates for the finite element method[END_REF]. Let h T be the diameter (see e.g. [START_REF] Scott | The Finite Element Method for Elliptic Problems[END_REF]) of the cell T and h E be the diameter of the facet E. The explicit residual estimator [START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF] on a cell T for the Poisson problems eqs. (2.2) and (2.4) is defined as

η 2 res,T := h 2 T f T + ∆u k 2 T + E∈E I ∩∂T 1 2 h E ∂ n u k E 2 E + E∈E N ∩∂T h E g E -∂ n u k 2 E , (A.1)
where f h and g h are the L 2 projections of f and g on V k respectively. In order to take into account inhomogeneous Dirichlet boundary conditions, we define in addition the Dirichlet oscillations. If E := Γ D ∩ T = ∅, then

osc 2 D,E := h E ∇ Γ (g E -u k ) 2 L 2 (E) , (A.2)
where ∇ Γ is the surface gradient and g E := π + T (g) is the L 2 projection of g onto V k+1

T [START_REF] Aurada | Each H 1/2 -stable projection yields convergence and quasi-optimality of adaptive FEM with inhomogeneous Dirichlet data in R d[END_REF]. The global residual estimator reads Abstract: This manuscript is concerned with a posteriori error estimation for the finite element discretization of standard and fractional partial differential equations as well as an application of fractional calculus to the modeling of the human meniscus by poroelasticity equations. In the introduction, we give an overview of the literature of a posteriori error estimation in finite element methods and of adaptive refinement methods. We emphasize the state-of-the-art of the Bank-Weiser a posteriori error estimation method and of the adaptive refinement methods convergence results. Then, we move to fractional partial differential equations. We give some of the most common discretization methods of fractional Laplacian operator based equations. We review some results of a priori error estimation for the finite element discretization of these equations and give the state-of-the-art of a posteriori error estimation. Finally, we review the literature on the use of the Caputo's fractional derivative in applications, focusing on anomalous diffusion and poro-elasticity applications. The rest of the manuscript is organized as follow. Chapter 1 is concerned with a proof of the reliability of the Bank-Weiser estimator for three-dimensional problems, extending a result from the literature. In Chapter 2 we present a numerical study of the Bank-Weiser estimator, provide a novel implementation of the estimator in the FEniCS finite element software and apply it to a variety of elliptic equations as well as goal-oriented error estimation. In Chapter 3 we derive a novel a posteriori estimator for the L2 error induced by the finite element discretization of fractional Laplacian operator based equations.

In Chapter 4 we present new theoretical results on the convergence of a rational approximation method with consequences on the approximation of fractional norms as well as a priori error estimation results for the finite element discretization of fractional equations. Finally, in Chapter 5 we provide an application of fractional calculus to the study of the human meniscus via poro-elasticity equations.

Université Bourgogne Franche-Comté 32, avenue de l'Observatoire 25000 Besançon

1 . 1 7 2

 117 Introduction xiii A posteriori error estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv Adaptive refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxii Fractional Laplacians . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii Fractional diffusion in poro-elasticity . . . . . . . . . . . . . . . . . . . . . . . . . xxix Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxi 1 Reliability of the Bank-Weiser estimator in dimension three 1 Model problem and finite element discretization . . . . . . . . . . . . . . . . 2 1.2 A posteriori error estimators . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Implementation of the Bank-Weiser estimator in FEniCSx 9 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . .

  Using quadrature rules in dimensions d = 2 and d = 3 give respectively T ψ F,T = 2 3 |T | and T ψ F,T = 3 10 |T |. In addition, by Poincaré's inequality and scaling arguments we have ∇ψ F,T T C h T |T |, and ψ F,T T C |T |. (1.16)

Figure 1 . 1 :

 11 Figure 1.1: On the left: The triangulation after four refinement steps of the adaptive algorithm driven by the Bank-Weiser estimator. On the right: The convergence curves for the estimator and exact error for each adaptive refinement algorithm.

  We consider the Poisson problem with mixed Dirichlet and Neumann boundary conditions. Let Γ = Γ D ∪ Γ N be a partition of the boundary. We apply a Dirichlet boundary condition on Γ D and a Neumann boundary condition on Γ N . Let f ∈ L 2 (Ω), u D ∈ H 1/2 (Γ D ) and g ∈ L 2 (Γ N ) be known data. We seek a function u:

Figure 2 . 1 :

 21 Figure 2.1: Overall process of the Bank-Weiser estimator algorithm.

  the Neumann boundary condition (traction) data g ∈ L 2 (Γ N ) d and displacement field u : Ω → R d . The stress and strain tensors are defined by σ := 2µε(u) -p Id, (2.40a) ε(u) := 1 2 ∇u + (∇u) T . (2.40b) where Id is the d × d identity matrix and µ and λ are the Lamé coefficients. The weak form of this linear elasticity problem reads: find u in H 1 (Ω) of trace u D on Γ D and p ∈ L 2 (Ω) such that

  .43b) Similarly to eqs. (2.41a) and (2.41b) transposed to the discrete context, eqs. (2.43a) and (2.43b) have a unique solution. If we denote e := uu 2 and ε := p -p 1 the discretization error is measured by 2µ ∇e T + r T .

2. 7 . 1

 71 Indicative snippet of error estimation for Poisson equation using Bank-Weiser estimator We present here a snippet of DOLFIN Python code showing function to compute the error of a Poisson problem using the Bank-Weiser estimator. from dolfin import * import fenics_error_estimation def estimate(u_h): """Bank-Weiser error estimation procedure for the Poisson problem. Parameters -----------u_h: dolfin.Function Solution of Poisson problem.

  Higher order space element_f = FiniteElement("DG", triangle, 2) # Low order space element_g = FiniteElement("DG", triangle, 1) # Construct the Bank-Weiser interpolation operator according to the # definition of the high and low order spaces. N = fenics_error_estimation.create_interpolation(element_f, element_g) V_f = FunctionSpace(mesh, element_f) e = TrialFunction(V_f) v = TestFunction(V_f) f = Constant(0.0) # Homogeneous zero Dirichlet boundary conditions bcs = DirichletBC(V_f, Constant(0.0), "on_boundary", "geometric") # Define the local Bank-Weiser problem on the full higher order space n = FacetNormal(mesh) a_e = inner(grad(e), grad(v))*dx # Residual L_e = inner(f + div(grad(u_h)), v)*dx + \ inner(jump(grad(u_h), -n), avg(v))*dS # Local solves on the implied Bank-Weiser space. The solution is returned # on the full space. e_h = fenics_error_estimation.estimate(a_e, L_e, N, bcs) # Estimate of global error error = norm(e_h, "H10") # Computation of local error indicator. V_e = FunctionSpace(mesh, "DG", 0) v = TestFunction(V_e) eta_h = Function(V_e, name="eta_h") # By testing against v in DG_0 this effectively computes # the estimator on each cell. eta = assemble(inner(inner(grad(e_h), grad(e_h)), v)*dx) eta_h.vector()[:] = eta return eta_h

Figure 2 . 2 :

 22 Figure 2.2: L-shaped Poisson problem with linear elements: On top left the initial mesh used to start all the adaptive strategies. From top middle to bottom right, the adaptive meshes obtained after seven iterations of refinement strategies steered respectively by η res , η b bw , η zz , η 2,1 bw and η 4,2 bw

Figure 2 . 3 : 84 ∅Figure 2 . 4 :

 238424 Figure 2.3: L-shaped Poisson problem with linear elements: efficiencies of η k + ,k - bw and other estimators on the last mesh of an adaptively refined hierarchy.

Figure 2 . 9 :

 29 Figure 2.9: L-shaped goal-oriented Poisson problem with linear elements: On top left the initial mesh used to start all the adaptive strategies. From top middle to bottom right, the adaptive meshes obtained after seven iterations of refinement strategies steered by weighted estimators derived respectively from η 4,2 bw , η res , η b bw and η zz for both primal and dual problems.

Figure 2 .

 2 Figure2.10: L-shaped goal-oriented Poisson problem with linear elements: plot comparing convergence of some goal-oriented adaptive strategies driven by four different estimators. Expected rates for primal and dual problems (-0.5) and goal functional (-1) shown by triangle markers. Comparison with an indicative line representing the higher order approximation of the errors of each strategy and obtained using least squares method.

Figure 2 .

 2 Figure 2.11: L-shaped goal-oriented Poisson problem with linear elements: efficiencies of the dual weighted estimators derived from η k + ,k - bw

2. 8 . 1

 81 Indicative snippet of error estimation for linear elasticity equations using Poisson estimatorWe give here a snippet of DOLFIN Python code showing function to compute the error of a two-dimensional linear elasticity problem (discretized with Taylor-Hood element) using the Poisson estimator, based on our implementation of the Bank-Weiser estimator. import scipy.linalg as sp.linalg from dolfin import * import fenics_error_estimation def estimate(w_h, mu, lmbda): """ Parameters -----------w_h: dolfin.Function Solution of the linear elasticity problem. mu: float First Lamé coefficient. lmbda: float Second Lamé coefficient. Returns -------The error estimate on each cell of the mesh. """ mesh = w_h.function_space().mesh() u_h = w_h.sub(0) p_h = w_h.sub(1) # Vectorial high order space. X_element_f = VectorElement('DG', triangle, 3) # Scalar high order and low order spaces. S_element_f = FiniteElement('DG', triangle, 3) S_element_g = FiniteElement('DG', triangle, 2) # Construct the scalar projection matrix according to the definition # of the high and low order spaces. N_S = create_interpolation(S_element_f, S_element_g) # Construct the vectorial projection matrix as a block diagonal, each # block corresponding to a scalar problem. N_X = sp.linalg.block_diag(N_S, N_S) f = Constant((0., 0.)) X_f = FunctionSpace(mesh, X_element_f) e_X = TrialFunction(X_f) v_X = TestFunction(X_f) # Homogeneous zero Dirichlet boundary conditions. bcs = DirichletBC(X_f, Constant((0., 0.)), 'on_boundary', 'geometric') # Cell residual. R_T = f + div(2.*mu*sym(grad(u_h))) -grad(p_h) # Facet residual. n = FacetNormal(mesh) R_E = (1./2.)*jump(p_h*Identity(2) -2.*mu*sym(grad(u_h)), -n) # Local Poisson problem. a_X_e = 2.*mu*inner(grad(e_X), grad(v_X))*dx L_X_e = inner(R_K, v_X)*dx -inner(R_E, avg(v_X))*dS # Solve Poisson equation locally on implicit Bank--Weiser space. e_h = fenics_error_estimation.estimate(a_X_e, L_X_e, N_X, bcs) # Cell residual. rho_d = 1./(lmbda**(-1)+(2.*mu)**(-1)) r_T = rho_d*(div(u_h) + lmbda**(-1)*p_h) # Computation of local error indicator. V_e = FunctionSpace(mesh, 'DG', 0) v = TestFunction(V_e) eta_h = Function(V_e) # By testing against v in DG_0 this effectively computes the estimator # on each cell. eta = assemble(2.*mu*inner(inner(grad(e_h), grad(e_h)), v)*dx + \ rho_d**(-1)*inner(inner(eps_h, eps_h), v)*dx) eta_h.vector()[:] = eta return eta_h 2.8.2 Nearly-incompressible elasticity We consider the linear elasticity problem from [89] on the centered unit square domain Ω with homogeneous Dirichlet boundary conditions on Γ D = Γ (u D = 0). The first Lamé coefficient is set to µ = 100 and the Poisson ratio to ν = 0.3 and ν = 0.499. The problem data f is given by f

1 ∅Figure 2 . 12 : 13 ∅Figure 2 . 13 :

 121213213 Figure 2.12: Nearly-incompressible elasticity (ν = 0.3) problem with Taylor-Hood elements: efficiencies of the Poisson estimators derived from η k + ,k - bw and other estimators on the last mesh of an adaptively refined hierarchy.

Figure 2 . 14 :uFigure 2 . 15 :

 214215 Figure 2.14: Femur bone linear elasticity problem with Taylor-Hood elements: on the top, the three different regions of the boundary corresponding to different boundary conditions: the left dark gray region is the non-zero Neumann boundary, the middle light gray region is the zero Neumann boundary and the right dark gray region is the Dirichlet boundary. In the middle, the initial mesh. On the bottom, the last mesh after several steps of adaptive refinement.

Figure 2 . 16 :

 216 Figure 2.16: Strong scaling study on the University of Luxembourg Aion HPC. Wall times for primal linear system assembly, primal linear system solve and error estimation of a threedimensional Poisson finite element problem on the unit square, discretized with quadratic elements. 1/1 triangle represents ideal strong scaling.

  -u κ,p ) |T , ∀T ∈ T . (3.23) Finally, we can estimate the L 2 error on the cell T by taking the norm of the function e bw κ,T e bw κ,T L 2 (T ) ≈ u κ -u κ,p L 2 (T ) . (3.24) These heuristics are summarized in fig. 3.1.

Figure 3 . 1 :

 31 Figure 3.1: Summary of the computation of the fractional solution approximation and of the fractional Bank-Weiser solution.

Figure 3 . 3 :

 33 Figure 3.3: Two-dimensional product of sines test case: the Bank-Weiser estimator η bw κ in solid blue line is compared to the exact error in dashed light blue line for s = 0.3 and s = 0.7.

Figure 3 . 4 :

 34 Figure 3.4: Two-dimensional product of sines text case: variation of the exact parametric errors with respect to the index l ∈ -M, N for three different fractional powers.

Figure 3 . 5 :

 35 Figure 3.5: Three-dimensional product of sines test case: the Bank-Weiser estimator η bw κ in solid blue line is compared to the exact error in dashed light blue line for various fractional powers.

Figure 3 . 6 :

 36 Figure 3.6: Two-dimensional checkerboard test case: for each fractional power we compare the values of the Bank-Weiser estimator η bwκ when uniform refinement is performed (dashed light lines) and when adaptive refinement is performed (solid lines).

  ] -0.35 -0.55 -0.75 -0.95 -1.00 Est. (unif.) -0.35 -0.55 -0.76 -0.95 -1.00 Est. (adapt.) -0.65 -0.84 -0.93 -0.97 -1.01 Table

Figure 3 . 7 :

 37 Figure 3.7: Two-dimensional checkerboard test case: meshes obtained after 10 steps of adaptive refinement steered by the Bank-Weiser estimator for s = 0.1, s = 0.5 and s = 0.9 from left to right.

Figure 3 . 8 :

 38 Figure 3.8: Three-dimensional checkerboard test case: finite element solution and mesh after six steps of adaptive refinement when s = 0.5. The unit cube domain (0, 1) 3 is truncated by the three planes passing through the point (0.25, 0.25, 0.25) and orthogonal to the vectors (1, 0, 0), (0, 1, 0) and (0, 0, 1) respectively.

Figure 3 . 9 :

 39 Figure3.9: Three-dimensional checkerboard test case: for each fractional power we compare the values of the Bank-Weiser estimator η bw κ when uniform refinement is performed (dashed light lines) and when adaptive refinement is performed (solid lines).

(4. 35 )

 35 In consequence, using eqs. (4.32), (4.34) and (4.35) we have|g λ (x -iη)| + |g λ (x + iη)| |g λ (x -iπ/4)| + |g λ (y + iπ/4)|. (4.36) So the integrand of eq. (4.25) is bounded uniformly in η by |g λ (x -iπ/4)| + |g λ (x + iπ/4)|. Let us show that |g λ (x -iπ/4)| + |g λ (x + iπ/4)| is integrable. First, following eq. (4.33)

(4. 47 )

 47 So eq. (4.45) is true and eq. (4.44) follows.From eq. (4.44) we deduce an upper bound for the function g λ . For any y ∈ R, we have g λ (y) = e 2sy 1 + e 2y λ -1

. 81 )

 81 Now, if N > M using again the arguments of eq. (4.51) we haveB 3 = κ +∞ l=N +1 |g λ (-lκ)| κ +∞ l=N +1 e -2slκ e -2sN κ 1 2s , (4.82) and since N = M + k with k > 0, B 3 e -2s(M +k)κ 1 2s e -2sM κ e -2skκ 1 eqs. (4.80), (4.81) and (4.83) the term B 3 is bounded by

Figure 4 . 1 :

 41 Figure 4.1: Plots of the function λ → λ s λ -s -Q N s (λ) for a wide range of values of λ and various values of N . From left to right, the fractional powers are s = 0.3, s = 0.5 and s = 0.7 respectively.

Figure 5 . 1 :

 51 Figure 5.1: (a) Schematic representation of the knee joint; (b) The different regions of the meniscus: posterior horn, central body, anterior horn and the three directions: vertical, radial and circumferential; (c) schematic representation of confined compression test set up, (d) example of the variability of the porosity within meniscal samples from µCT scans; (e) results show that the central body of the meniscus exhibits higher permeability values .

  classical linear model of transient flow and deformation of a homogeneous fully saturated elastic porous medium depends on an appropriate coupling of the fluid pressure and solid stress. A change in applied stress produces a change in fluid pressure or fluid mass and a change in fluid pressure or fluid mass is responsible for a change in the volume of the porous material. The coupling term affects only the hydrostatic part of the stress tensor. The stress tensor can be written as follows: σ = 2Gε + λtrace(ε)I -αpI (5.11)

Figure 5 . 2 :

 52 Figure 5.2: (a) Confined compression test set up. Cylindrical specimens are extracted from a human meniscus. (b) The cylinder is then confined inside a glass chamber with a porous base and compressed by a metal bar. During compression, the water in the meniscal sample flows from the base of the cylinder.(b) Schematic representation of the boundary conditions used to solve the boundary value problem.

Figure 5 . 3 :

 53 Figure 5.3: Examples of fittings of relaxation test data (only one of the 5-steps test is reported) with eq. (5.10) for a medial human meniscus (TK11) -central body (CB). a) b) and d) samples extracted in radial direction and c) in vertical direction. Table2contains the parameters H A and permeability k of all the 5-steps.

Figure 5 . 4 :

 54 Figure 5.4: Evolution of the decrease in weight during the confined compression test which has been correlated with the amount of the fluid discharged upon compression of the tissue. The fittings relation is given in eq. (5.26). The graphs shown are related to samples extracted from a lateral meniscus (tk16). (a) Sample from the central body (B) and along the circumferential direction (C). (b) Sample from the posterior horn (P) and along the radial direction (R). (c) Sample from the central body (B) and along the vertical direction (V). (d) Sample from the posterior horn (P) and along the vertical direction (V). (e) Sample from the central body (B) and along the radial direction (R). (f) Sample from the anterior horn (A) and along the radial direction (R).

Figure 5 . 5 :

 55 Figure 5.5: Graphical explanation of the region (central body) /direction (vertical, radial and circumferential) corresponding to the parameters highlighted in table 5.5.

  λ β are higher for the central body of the meniscus with respect of the anterior and posterior horns. Mean values of λ β , β in the central body are λ β = 5.5443 × 10 -10 m 4 Ns 1-β , β = 0.0434, while in the posterior and anterior regions are λ β = 2.851 × 10 -10 m 4 Ns 1-β , β = 0.0326 and λ β = 1.2636 × 10 -10 m 4 Ns 1-β , β = 0.0232 respectively.

Figure 5 .

 5 Figure 5.6a shows the normalized flow rate ( jf ) calculated considering eq. (5.1) in the case of a constant value of gradient of pressure ∇p = 3 × 10 8 Pa m , λ β = 5.5443 × 10 -10 m 4 N× 1-β and the three values of β = 0.0434, 0.0326, 0.0232 obtained for the central body, posterior and anterior horns respectively. It can be noted that even a small value of the time fractional derivative β affects the fluid flow response. β = 0 is equivalent to considering the classical Darcy law, i.e. fluid flow rate is constant in time. From fig. 5.6ait is possible to observe that the higher the value of β the faster the decrease in the fluid flow rate is. Furthermore, we analyzed the evolution of the fluid flow rate in the central body, posterior and anterior horns using the mean values of the parameters β, λ β above. Figure5.6b shows the response of the normalized flux in the three regions. It can be noted that as the value of the anomalous permeability λ β in the central body is about four times higher than in the anterior horn, the value of flux is higher in the central body and lowest in the anterior region. Moreover, as the value of β is larger compared to the anterior and posterior horns, the fluid flux is faster (i.e. decrease in fluid flow rate) in the central body. It would be also important to consider how the fluid flow rate evolves in different directions, i.e. radial, vertical or circumferential. In this regard, given the paucity of data related to anterior and posterior horns, it is only possible to analyze the central body of the meniscus. Mean values of β, λ β in the circumferential, radial and vertical directions are β = 0.0421, λ β = 6.3924 × 10 -10 m 4 Ns (1-β) , β = 0.0488, λ β = 5.9196 × 10 -10 m 4 Ns (1-β) and β = 0.0393, λ β = 4.321 × 10 -10 m 4Ns (1-β) , respectively. Figure5.6c pictures the evolution

  Figure 5.6a shows the normalized flow rate ( jf ) calculated considering eq. (5.1) in the case of a constant value of gradient of pressure ∇p = 3 × 10 8 Pa m , λ β = 5.5443 × 10 -10 m 4 N× 1-β and the three values of β = 0.0434, 0.0326, 0.0232 obtained for the central body, posterior and anterior horns respectively. It can be noted that even a small value of the time fractional derivative β affects the fluid flow response. β = 0 is equivalent to considering the classical Darcy law, i.e. fluid flow rate is constant in time. From fig. 5.6ait is possible to observe that the higher the value of β the faster the decrease in the fluid flow rate is. Furthermore, we analyzed the evolution of the fluid flow rate in the central body, posterior and anterior horns using the mean values of the parameters β, λ β above. Figure5.6b shows the response of the normalized flux in the three regions. It can be noted that as the value of the anomalous permeability λ β in the central body is about four times higher than in the anterior horn, the value of flux is higher in the central body and lowest in the anterior region. Moreover, as the value of β is larger compared to the anterior and posterior horns, the fluid flux is faster (i.e. decrease in fluid flow rate) in the central body. It would be also important to consider how the fluid flow rate evolves in different directions, i.e. radial, vertical or circumferential. In this regard, given the paucity of data related to anterior and posterior horns, it is only possible to analyze the central body of the meniscus. Mean values of β, λ β in the circumferential, radial and vertical directions are β = 0.0421, λ β = 6.3924 × 10 -10 m 4 Ns (1-β) , β = 0.0488, λ β = 5.9196 × 10 -10 m 4 Ns (1-β) and β = 0.0393, λ β = 4.321 × 10 -10 m 4Ns (1-β) , respectively. Figure5.6c pictures the evolution

Figure 5 .

 5 Figure 5.6: a. Normalized fluid flow rate ( jf ) calculated considering eq. (5.1) in the case of a constant value of gradient of pressure ∇p = 3 × 10 8 P a/m, λ β = 5.5443 × 10 -10 m 4 N s 1-β and three values of β = 0.0434, 0.0326, 0.0232 obtained for the central body, anterior and posterior horns. b. Normalized fluid flow rate ( jf ) for the three regions (central body, anterior and posterior horns) calculated considering eq. (5.1) in the case of a constant value of gradient of pressure ∇p = 3 × 10 8 P a/m and with λ β = 5.5443 × 10 -10 m 4 N s 1-β , β = 0.0434 in the central body, λ β = 2.851 × 10 -10 m 4 N s 1-β , β = 0.0326 in the posterior horn and λ β = 1.2636 × 10 -10 m 4 N s 1-β , β = 0.0232 in the anterior horn. c. Evolution of the normalized fluid flow rate in circumferential, radial and vertical directions of the central body portion of the meniscus with β = 0.0421, λ β = 6.3924 × 10 -10 m 4 N s 1-β for the circumferential direction, β = 0.0488, λ β = 5.9196 × 10 -10 m 4 N s 1-β for the radial direction and β = 0.0393, λ β = 4.321 × 10 -10 m 4N s 1-β for the vertical direction.

Figure 5 . 7 :

 57 Figure 5.7: Evolution of the adimensional pore pressure (p = p(z, t)/P A in eq. (5.21)) throughout the length of the sample in the z axis for different values of β = 0, 0.1, 0.5 at a time of 15s.

A. 2 2 T + ρ d r T 2 T

 222 Linear elasticity equationsThe residual estimator for the linear elasticity problem eqs. (2.41a), (2.41b), (2.43a) and (2.43b) is given by η 2 res,T := ρ T R T R T , r T and R E are respectively defined in eqs. (2.44a) to (2.44c) and the constants ρ T , ρ d and ρ E are given byρ T := h T (2µ) -1/2 2 , ρ d := λ -1 + (2µ) -1 -1 , ρ E := h E (2µ)Estimation d'erreur a posteriori pour l'approximation de problèmes Laplaciens fractionnaires et applications en poro-élasticité Mots clés: Méthodes éléments finis, Estimation d'erreur a posteriori, Equations aux dérivées partielles fractionnaires, Poroélasticité Résumé: Ce manuscrit traite d'estimation d'erreur a posteriori pour la discrétisation d'équations aux dérivées partielles standard et fractionnaires par les méthodes éléments finis ainsi que de l'application de l'analyse fractionnaire à la modélisation du ménisque humain par les équations de poro-élasticité. Dans l'introduction, nous donnons un aperçu de la littérature sur l'estimation d'erreur a posteriori pour les méthodes éléments finis et des méthodes de raffinement adaptatif. Nous insistons particulièrement sur l'état de l'art de la méthode d'estimation d'erreur a posteriori de Bank-Weiser et sur les résultats de convergence des méthodes adaptatives. Ensuite, nous nous intéressons aux équations aux dérivées partielles fractionnaires. Nous présentons certaines méthodes de discrétisation d'équations basées sur l'opérateur Laplacien fractionnaire et donnons l'état de l'art sur l'estimation d'erreur a posteriori. Finalement, nous donnons un aperçu de la littérature concernant les applications de la dérivée fractionnaire au sens de Caputo en nous concentrant sur le phénomène de diffusion anormale et les applications en poro-élasticité. Title: A posteriori error estimation for finite element approximations of fractional Laplacian problems and applications to poro-elasticity Keywords: Finite element methods, A posteriori error estimation, Fractional partial differential equations, Poroelasticity

  Primal problem statement and finite element discretization . . . . . . . . . . 2.3 The Bank-Weiser estimator . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 The global error equation . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 The local Bank-Weiser space and the Bank-Weiser estimator . . . . 2.3.3 A particular example . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Algorithms and implementation details . . . . . . . . . . . . . . . . . . . . . 2.4.1 Method outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Computational details . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.3 Additional remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Adaptive mesh refinement . . . . . . . . . . . . . . . . . . . . . . . iii First of all, I would like to acknowledge the support of the ASSIST research project of the University of Luxembourg. This work has been prepared in the framework of the DRIVEN TWINNING project funded by the European Union's Horizon 2020 Research and Innovation programme under Grant Agreement No. 811099. I would like to thank Prof. Andrea Bonito and Prof. Martin Vohralík to have accepted to be my dissertation examiners and to be part of my defence jury. It is a great honor for me to share my research with them.

. . . . . 2.1.2 Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.4 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2

  .14) Note, due to eq. (2.13), the matrix of S in the couple of basis (B + T , B + T ) is the identity matrix of size d + × d + . Consequently, if we denote G the matrix of L T in the basis B + T and G the matrix of L T in the basis B -

T

, we have

  .6 are very different from fig. 2.4, many Bank-Weiser estimators are now underestimating the error. The most efficient estimator is η 2,1 bw closely followed by the bubble Bank-Weiser estimator η b bw . As

		k -	k +	1	2	3	4
			0	0.83 1.06 1.08 1.14
			1	∅ 0.94 1.21 1.34	η res	2.84
			2	∅	∅	0.0 0.55	η b bw	1.24
			3	∅	∅	∅ 0.23	η zz	0.91
	Figure 2.5: Mixed boundary conditions L-shaped Poisson problem with linear elements:
	efficiencies of η	k + ,k -bw	and other estimators on the last mesh of an adaptively refined hier-
	archy.
		k -	k +	1	2	3	4
			0	0.57 0.86 0.83 0.97
			1	∅ 1.05 1.37 1.43	η res	5.91
			2	∅	∅ 0.62 0.78	η b bw	1.17
			3	∅	∅	∅	0.3	η zz

∅ Figure 2.6: Mixed boundary conditions L-shaped Poisson problem with quadratic elements: efficiencies of η k + ,k - bw

Table 3 . 3 :

 33 Two-dimensional product of sines test case: efficiency indices of the Bank-Weiser estimator for various fractional powers.

  biphasic model depends essentially on two parameters: the aggregate modulus H A and the averaged axial permeability λ 0 . The linear fractional poroelasticity model of Biot

	Notation Name	Expression	Unit
	K	drained bulk modulus		Pa
	G	drained shear modulus		Pa
	K u	undrained bulk modulus		Pa
	α	Biot's coefficient		Unitless
	B	Skempton's coefficient	(K u -K)/(αK u ) Unitless
	k	averaged intrinsic permeability		m 2
	µ	fluid viscosity		N • s • m -2
		averaged axial permeability
	λ 0	-compressible case

Table 5 .

 5 1: Parameters of the linear (fractional) Biot model.

	Notation Name	Expression	Unit
	γ	ratio solid volume/fluid volume		Unitless
	K s	solid phase bulk modulus		Pa
	G s	solid phase shear modulus		Pa
	k	intrinsic permeability		m 2
	µ	fluid viscosity		N • s • m -2
	λ 0	averaged axial permeability	µ/((1 + γ) 2 k) N • s • m -4
	H A	aggregate modulus	(4G s + 3K s )/3 Pa

Table 5 . 2

 52

: Parameters of the Biphasic model.

Table 5 .

 5 4: Complete set of the 5 step test results of the central portion of the medial meniscus shown in fig.5.3 for the three directions (vertical, radial, circumferential). Note that the expected trend of decreasing permeability with increasing strain was visible only for the vertical direction (TK11-CB-Vert). Root mean square error (RMSE) values are high in some cases showing that the model does not always fit the experimental behavior.

			2	0.180	0.858 × 10 -12 0.0003
		2	4	0.135	2.803 × 10 -12 0.0002
	TK11-CB-Circ	3	6	0.032	0.410 × 10 -12 0.0011
		4	8	0.122	0.953 × 10 -12 0.0004
		5	10	0.058	1.714 × 10 -12 0.0004
		1	2	-	-	-
		2	4	0.044	1.572 × 10 -12 0.0005
	TK11-CB-Rad	3	6	0.007	3.383 × 10 -12 0.1533
		4	8	0.037	0.870 × 10 -12 0.0003
		5	10	0.030	0.525 × 10 -12 0.0008
		1	2	0.0173	1.0190 × 10 -12 0.0023
		2	4	0.0116	0.9400 × 10 -12 0.0100
	TK11-CB-Vert	3	6	0.0126	0.8100 × 10 -12 0.0120
		4	8	0.0230	0.4760 × 10 -12 0.0025
		5	10	0.0173	0.4850 × 10 -12 0.0080

Table 5 .

 5

5: Results of the fittings. First column shows the name of the sample, second column β the order of derivative, third column λ β is the anomalous permeability, forth column specifies if the sample is extracted from a medial (M) or lateral (L) meniscus.

Merci à mes parents, je sais combien de sacrifices vous avez fait pour que j'ai la chance de faire ce qui me plaît aujourd'hui. Je profite de cette occasion pour vous dire à quel point cela a été important pour moi.Merci à mon frère et ma soeur pour avoir été présents, chaque moment passé avec vous deux m'est précieux.Pour finir, je voudrais remercier celle qui a accepté de me suivre dans cette aventure et dans celles qui vont suivre, celle qui a partagé mes joies et mes peines pendant toutes

The STL model of the femur bone can be found at https://3dprint.nih.gov/discover/3dpx-000168 under a Public Domain license.

Acknowledgments APPENDIX B. THE ZIENKIEWICZ-ZHU ESTIMATOR

Bibliography

Supplementary material

The two versions of the code can be found at • https://github.com/jhale/fenicsx-error-estimation (DOLFINx version),

• https://github.com/rbulle/fenics-error-estimation (DOLFIN version).

A simplified version of the code (LGPLv3) used to produce the results in this paper is archived at https://doi.org/10.6084/m9.figshare.10732421. A Docker image [START_REF] Hale | Containers for Portable, Productive, and Performant Scientific Computing[END_REF] is provided in which this code can be executed.

Choose a tolerance ε > 0 Choose an initial mesh T n=0 Choose κ such that ε s (κ) f L 2 < ε Generate the rational approximation Q computation of the local Bank-Weiser solutions on the cells are independent from each other. Our error estimation strategy combines these advantages and is fully parallelizable both with respect to the parametric problems and local estimators computation. An example of error estimation and adaptive refinement algorithm based on our method is shown in fig. 3.2.

The algorithm presented in fig. 3.2 is based on three loops: one While loop and two For loops. The While loop is due to the adaptive refinement procedure and can not be parallelized. However, the two For loops are fully parallelizable and this parallelization can be highly advantageous for large three-dimensional problems.

Note that there is no guarantee that the mesh we obtain at the end of the main While loop in fig. 3.2 is optimal for all the parametric problems. For some of the parametric solutions without boundary layers the mesh is certainly over-refined. An alternative approach could be to compute the L 2 norms of the parametric Bank-Weiser solutions e bw l,T in order to derive parametric Bank-Weiser estimators and refine the meshes independently for each parametric problem. This would require the storage of a possibly different mesh for each Chapter 5 Application: poro-elasticity of the human meniscus

Abstract

In this study, we observe that the poromechanical parameters in human meniscus vary spatially throughout the tissue. The response is anisotropic and the porosity is functionally graded. To draw these conclusions, we measured the anisotropic permeability and the "aggregate modulus" of the tissue, i.e. the stiffness of the material at equilibrium, after the interstitial fluid has ceased flowing. We estimated those parameters within the central portion of the meniscus in three directions (i.e., vertical, radial and circumferential) by fitting an enhanced model on stress relation confined compression tests. In fact, we noticed that a classical biphasic model was not sufficient to reproduce the observed experimental behavior. We propose to adopt a poroelastic model based on the assumption that the fluid flow inside the human meniscus is described by a fractional porous medium equation analogous to Darcy's law which involves fractional operators. The fluid flux is then time-dependent for a constant applied pressure gradient (in contrast with the classical Darcy's law which describes a time independent fluid flux relation). We show that a fractional poroelastic model is well-suited to describe the flow within the meniscus and identify the associated parameters (i.e., the order of the time derivative and the permeability). The Zienkiewicz-Zhu estimator

The Zienkiewicz-Zhu estimator is a gradient recovery estimator based on an averaging technique introduced in [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF]. This estimator belongs to a general class of recovery estimators, see [START_REF] Cai | Recovery-Based Error Estimator for Interface Problems: Conforming Linear Elements[END_REF][START_REF] Cai | Improved ZZ a posteriori error estimators for diffusion problems: Conforming linear elements[END_REF][START_REF] Zhang | Recovery type a posteriori error estimates in finite element methods[END_REF] for recent surveys and a reformulation of the recovery procedure in an H(div)-conforming space that has superior performance for problems with sharp interfaces. Despite the fact that some recovery estimators, especially when based on least squares fitting, are available for higher order finite elements (see for example [START_REF] Zienkiewicz | The superconvergent patch recovery (SPR) and adaptive finite element refinement[END_REF]) we only consider the original estimator, defined for a piecewise linear finite element framework. Given the finite element solution u 1 ∈ V 1 the numerical flux ρ 1 := ∇u 1 is a piecewise constant vector field. For each vertex χ ∈ N in the mesh we denote ω χ the domain covered by the union of cells T having common vertex χ. The code in the supplementary material contains a prototype implementation of the Zienkiewicz-Zhu estimator in FEniCS. We have implemented the local recovered flux calculation in Python rather than C++, so the runtime performance is far from optimal.