You are an expert at summarizing long articles. Proceed to summarize the following text: 
flow inside a porous medium finds many applications in natural and engineering systems . subsurface flows  @xcite , erosion  @xcite , fuel cells  @xcite , and filtration systems  @xcite are a few examples of physical processes that are governed by low reynolds number porous media flow . accurately resolving such flows is essential for modelling transport  @xcite , mixing  @xcite , anomalous diffusion  @xcite , breakthrough curves  @xcite , carbon dioxide sequestration  @xcite , uncertainty quantification  @xcite , chemical reactions  @xcite , and many other applications . in most of these works , in particular , those that focus on numerics and modelling , the porous medium is assumed to be two dimensional . therefore , there is a strong need for efficient and accurate numerical methods to simulate complex two - dimensional porous media flow . in this paper , we focus on solving the two - dimensional incompressible stokes equations in the geometry illustrated in . the incompressible stokes equations have been solved in porous domains resembling using finite differences  @xcite , finite volumes  @xcite , algebraic multigrid  @xcite , lattice - boltzmann  @xcite , and continuous time random walks  @xcite . however , each of these methods struggle with either enforcing incompressibility , meshing the geometry , or obtaining high - order accuracy . boundary integral equations ( bies ) , on the other hand , are a powerful method that addresses each of these challenges . in particular , bies are naturally adaptive , allow for high - order approximations , automatically enforce the incompressibility constraint , and enable the computation of the velocity field without computing the pressure . another advantage of a bie is that the unknown two - dimensional velocity field is written in terms of a one - dimensional integral over the boundary of the domain . this dimension reduction implies that only the boundary , which is a collection of one - dimensional closed curves , needs to be meshed . furthermore , a natural extension of this work is to solve the three - dimensional stokes equations while only discretizing the two - dimensional boundary of the geometry . other groups have used integral equations to simulate porous media by using a darcy approximation  @xcite , or by coupling a penalization method with a volume integral equation  @xcite . in contrast , we will use a boundary integral equation to solve the incompressible stokes equations . the main disadvantage of bies is the need to solve a dense linear system of equations . a direct solver based on gaussian elimination requires @xmath3 operations  where @xmath4 denotes the number of unknowns  which is impractical for realistic geometries . it also requires the assembly and storage of the dense matrix , leading to an @xmath5 memory requirement . iterative solvers such as gmres , on the other hand , only rely on ( accelerated ) matrix - vector multiplications and have an @xmath6 memory requirement , but can suffer from a poor convergence behavior due to the ill - conditioning of the linear system . in the recent decade , there have been two main directions of research that have lead to faster algorithms for solving discretized bies . a first approach focuses on developing integral equations that result in linear systems that are both well - conditioned and require a mesh - independent number of gmres iterations to reach convergence  @xcite . these properties generally follow from second - kind integral equation methods with compact integral operators . because of the mesh - independent number of gmres iterations , the asymptotic cpu time scales with the cost of a single dense matrix - vector multiplication . using fast summation methods such as the fast multipole method ( fmm )  @xcite , ewald summation  @xcite , or tree - code methods  @xcite , the cost of a matrix - vector multiplication can be reduced from @xmath5 to either @xmath7 or @xmath6 . while this approach can result in an optimal asymptotic complexity , complexities of the geometry , such as nearly touching pores and regions of high curvature , often still results in a large number of gmres iterations  @xcite . in addition , the double - layer potential formulation , which requires a mesh - independent number of gmres iterations , has less regularity than the single - layer potential  @xcite , which requires a mesh - dependent number of gmres iterations . to reduce the number of iterations , a preconditioner can be incorporated . while simple techniques such as block - diagonal preconditioners are easy to construct , they are not always capable of reducing the number of iterations significantly . an incomplete list of more sophisticated preconditioners , which vary in construction time and efficacy , include additive or multiplicative schwarz  @xcite , bpx  @xcite , sparse approximate inverses  @xcite , approximate lu decompositions  @xcite , clustering  @xcite , integral equations of opposite order  @xcite , calderon identities  @xcite , or multigrid  @xcite . in the case of multigrid preconditioners , specialized smoothers must be developed since standard smoothers like jacobi or gauss - seidel will smooth the low , rather than the high , frequency components . examples of such smoothers may require developing approximate inverses of the integral operator  @xcite involving , for example , the laplace - beltrami operator  @xcite . a second approach targets the use of fast direct solvers that exploit the structure of the linear system that arises when discretizing bies . the structure can be loosely described as matrices with hierarchical off - diagonal blocks being numerically low - rank . these matrices are formalized as , with an increasing level of complexity , hierarchically off - diagonal low - rank ( hodlr ) matrices  @xcite , hierarchically semi - separable ( hss ) matrices  @xcite , @xmath8-matrices  @xcite , and @xmath0-matrices  @xcite . fast methods for efficiently constructing and storing approximate inverses of these matrices have been developed in recent years  @xcite , mainly for the aforementioned hodlr and hss matrices . moreover , once the approximate inverse , which only depends on the geometry , is computed and stored , it can be applied with a minimal amount of cpu time to multiple right - hand sides . therefore , direct solvers are extremely useful when applying a temporal discretization to a partial differential equation  @xcite , or when considering scattering for many different incident waves  @xcite . we blend the two aforementioned approaches in this paper : we present the inverse fast multipole method ( ifmm )  @xcite as an inexact fast direct solver based on @xmath0-matrices and apply it as a preconditioner in a fast multipole accelerated gmres solver  @xcite . the preconditioner has a tunable accuracy @xmath1 and its computational cost scales almost linearly with the problem size as @xmath9 . we will demonstrate that the ifmm is capable of significantly reducing the number of iterations and the cpu time . the outline of this paper is as follows . discusses the fluid model and layer potential formulation . the inverse fast multipole method ( ifmm ) is subsequently introduced in , focusing on the main features for achieving an efficient preconditioner . numerical benchmarks are presented in . we report results for both the unpreconditioned system , and for preconditioned systems where block - diagonal and ifmm preconditioners ( at several accuracies ) are used . summarizing conclusions are finally drawn in . we are interested in simulating an incompressible newtonian fluid in the geometry illustrated in . we consider scales where the dimensionless reynolds number is small so that the fluid is assumed to be stokesian . if we assumed a homogeneous porosity , the much simpler darcy equations could be solved . however , we do not make this assumption and we focus on accurately and efficiently solving the incompressible stokes equation with a no - slip boundary condition on each pore . we also assume that the flow is two - dimensional and that there are no external body forces such as gravity . we avoid the stokes paradox by bounding the porous region by the boundary @xmath10 . to resemble a pipe flow , we let @xmath10 be a mollification of the rectangle @xmath11 \times [ -h , h]$ ] ( see ) . by smoothing the corners of @xmath10 , specialized quadrature for integral equations is avoided . each inner boundary , denoted by @xmath12 , @xmath13 , is a circle of variable radii . the area enclosed by the outer boundary is denoted by @xmath14 , and the area enclosed by each pore by @xmath15 . therefore , the geometry is given by @xmath16 and its boundary is @xmath17 .    with this setup , by defining the fluid velocity @xmath18 , the pressure @xmath19 , the viscosity @xmath20 , and a prescribed flow on the boundary @xmath21 , the governing equations are @xmath22 we impose boundary conditions that correspond to a pipe flow at the intake and outtake , and no slip on the boundary of the pores @xmath23 where @xmath24 sets the velocity scale . furthermore , we assume that @xmath25 from this point onwards . a bie formulation of has several advantages over the differential form . the pressure need not be computed , the incompressibility constraint is automatically satisfied , and the velocity is written in terms of a layer potential involving an unknown density function defined only on @xmath26 . two popular choices for the layer potential are the single - layer and double - layer potentials @xmath27({\mathbf{x } } ) & = \frac{1}{4\pi } \int_{\gamma } \left (     -\log \rho \mathbf{i } + \frac{{\mathbf{r}}\otimes { \mathbf{r}}}{\rho^{2}}\right )     { \boldsymbol{\sigma}}({\mathbf{y } } ) ds_{{\mathbf{y } } } , \qquad & & { \mathbf{x}}\in \omega , \\    \mathcal{d}[{\boldsymbol{\sigma}}]({\mathbf{x } } ) & = \frac{1}{\pi}\int_{\gamma }    \frac{{\mathbf{r}}\cdot { \mathbf{n}}}{\rho^{2 } } \frac{{\mathbf{r}}\otimes { \mathbf{r}}}{\rho^{2 } }    { \boldsymbol{\sigma}}({\mathbf{y } } ) ds_{{\mathbf{y } } } , & & { \mathbf{x}}\in \omega,\end{aligned}\ ] ] respectively , where @xmath28 , @xmath29 , @xmath30 is the unit outward normal , and @xmath31 is a density function defined on @xmath26 . given our dirichlet boundary condition , the double - layer potential formulation results in a second - kind integral equation whose condition number is guaranteed to be mesh - independent . however , the single - layer potential requires less smoothness of @xmath31 since  @xcite @xmath32 are continuous . moreover , this result holds even under the relaxed assumption that @xmath33 is a bounded lipschitz domain . another advantage of the single - layer potential is that it has full - rank , as opposed to the double - layer potential which has a rank three null space  @xcite for each connected component of @xmath26 . while this rank deficiency can be removed , it is unclear if this is compatible with the ifmm . moreover , for complex geometries , the double - layer potential may still require many gmres iterations  @xcite , but in this case , we expect that applying the ifmm preconditioner will result in a computational speedup . given the additional regularity and full - rank of the single - layer potential , we choose to represent the velocity @xmath34 as @xmath35({\mathbf{x } } ) = \frac{1}{4\pi } \int_{\gamma } \left (     -\log \rho \mathbf{i } + \frac{{\mathbf{r}}\otimes { \mathbf{r}}}{\rho^{2}}\right )     { \boldsymbol{\sigma}}({\mathbf{y } } ) ds_{{\mathbf{y } } } , \quad { \mathbf{x}}\in \omega ,    \label{eq : stokesslp}\end{aligned}\ ] ] where @xmath31 is an unknown density function . taking the limiting value of as @xmath36 tends to @xmath26 , the density function @xmath31 must satisfy the first - kind integral equation  @xcite @xmath37 it is known that upon discretizing , the result is an ill - conditioned linear system that requires a mesh - dependent number of gmres iterations . however , the linear system is amenable to the solvers for structured matrices . once the unique solution of is computed , the density function @xmath31 is substituted into to evaluate the velocity @xmath18 at any point @xmath38 . we adopt a collocation method to discretize . first , each interior curve @xmath12 , @xmath13 is parameterized and discretized at @xmath39 points , and the outer wall is discretized at @xmath40 points . the resulting discretization points are @xmath41 , @xmath42 . we enforce the bie at these @xmath43 collocation points by requiring @xmath44 then , quadrature is applied , where the quadrature nodes coincide with the discretization points @xmath41 . this results in the dense linear system @xmath45 where @xmath46 is the kernel of the integral operator in . the weights @xmath47 are the product of the arclength term with appropriate quadrature weights . we integrate the weak logarithmic singularity by applying a sixth - order quadrature rule of kapur and rokhlin  @xcite . this quadrature rule is identical to the trapezoid rule , except that the weights are modified at the six nodes to the left and right of the singularity . therefore , the quadrature nodes coincide with the discretization points @xmath41 , @xmath48 . since there are only @xmath49 quadrature weights that are modified for each @xmath41 , these quadrature rules are compatible with the fmm . in addition , in contrast to alpert quadrature rules  @xcite , additional quadrature nodes are not introduced , so no interpolation is required . given this quadrature formula ,   results in a dense @xmath50 linear system that , for brevity , we write as @xmath51 the matrix @xmath52 is a structured matrix with numerical low - rank off - diagonal blocks . gmres can be used to iteratively solve  , and each matrix - vector multiplication can be done with @xmath6 complexity by using the fmm . however , the number of required gmres iterations will be large , and depend on the mesh size . therefore , we will use the structure of @xmath52 to develop an efficient preconditioner for gmres . upon computing the density function @xmath31 at the points @xmath53 , the velocity field needs to be computed at points @xmath38 . this can be done with the trapezoid rule @xmath54 where @xmath55 is an arclength spacing of @xmath26 at @xmath56 . the trapezoid rule guarantees spectral accuracy , but this asymptotic convergence rate is delayed when @xmath36 is too close to @xmath26 . ] . there are methods to guarantee a uniform error for all @xmath57  @xcite , but this is not the focus of this work . instead , we focus on efficiently solving for the density function @xmath31 . the inverse fast multipole method ( ifmm )  @xcite is an inexact fast direct solver for @xmath0-matrices that can be applied as a preconditioner for gmres . it has successfully been used to precondition matrices arising from integral operators  @xcite and radial basis interpolation  @xcite . the main features of the algorithm are concisely described in the following subsections ; the reader is referred to @xcite for a more detailed description . although a 2d problem is considered in this paper , the ifmm is applicable to 3d problems as well ( see  @xcite for examples ) . a variety of fast direct solvers for dense linear systems such as have been developed in recent years . these methods are based on a multilevel , hierarchical decomposition of the physical domain of interest , from the root down to the leaf level  @xmath58 ( e.g. , by means of a quadtree in 2d or an octree in 3d , either uniform or adaptive ) . this is subsequently exploited to represent the matrix @xmath52 in a compressed format . the simplest of these methods assume that all off - diagonal blocks are low - rank ( only self - interactions are considered to be of full - rank ) , leading to so - called hierarchically off - diagonal low - rank ( hodlr )  @xcite and hierarchically semi - separable ( hss ) matrices  @xcite . in the latter , a nested approach is used ( i.e. ,  the low - rank basis at a certain hierarchical level is constructed using the low - rank basis at the child level ) , while this is not the case in the former . given this assumption on the matrix structure , the inversion can be performed exactly . direct @xmath6 algorithms for hodlr matrices can be found in  @xcite , while fast solvers for hss matrices have been presented by , among others , martinsson and rokhlin  @xcite , chandrasekaran et al . @xcite , xia et al . @xcite , and gillman et al .  @xcite . the assumption that all off - diagonal blocks are low - rank is a major drawback for the aforementioned hodlr and hss matrices . indeed , the rank is in that case likely to grow in an unbounded fashion for 2d and 3d problems with a complex geometry  @xcite , hence jeopardizing the computational efficiency of the algorithms . this drawback can be circumvented by adopting the more general framework of @xmath8- or @xmath0-matrices ( non - nested vs.  nested ) , in which only non - neighboring interactions are replaced by low - rank approximations . this is the approach followed in the ifmm  @xcite . it is worth noting that the algorithm presented in  @xcite has recently been extended by corona et al . @xcite to obtain a direct solver for hss matrices that scales for 2d problems as @xmath6 . we note that since our problem is two - dimensional , a hss approximation of would potentially have bounded rank . however , to extend the work to three dimensions , we choose to use an @xmath0 framework . an @xmath0-matrix @xmath59 is weak hierarchical matrix ( in which only non - neighboring , well - separated interactions are assumed to be low - rank ) with a nested basis that approximates the original matrix @xmath52 . it can be expressed as  @xcite @xmath60 @xmath59 is thus decomposed into a block sparse matrix , @xmath61 , containing all the neighboring and self - interactions at the leaf level @xmath58 , and a low - rank term , @xmath62 , that characterizes all well - separated interactions . @xmath63 and @xmath64 are interpolation and anterpolation operators , respectively , while @xmath65 is an @xmath0-matrix of a smaller size . is applied recursively to @xmath65 until the top level . ] , where there are no more well - separated interactions , is reached . the ifmm is a novel algorithm for the approximate inversion of an @xmath0-matrix , and consists of two crucial steps : extended sparsification and compression of fill - ins that appear throughout the elimination . the hierarchical structure of @xmath59 can be exploited to represent @xmath66 as an extended sparse system , which is more attractive from a computational point of view . introducing auxiliary variables @xmath67 and @xmath68 ( corresponding to multipole and local coefficients in the fmm ) leads to the following extended system : @xmath69 in , the sparse blocks @xmath61 , @xmath63 , and @xmath70 are situated in the top left corner of the matrix , while the remaining dense block @xmath65 is pushed to the bottom right corner of the matrix . as the latter is an @xmath0-matrix as well , the same technique can be applied recursively to further extend and sparsify the system . one finally obtains the following system of equations , which is equivalent to @xmath66 : @xmath71 the variables and equations are ordered in such a way that a sparse matrix with a symmetric fill - in pattern is obtained . note that this extended system is only moderately larger than the original system , since the dimensions of the auxiliary variables ( characterizing the low - rank interactions ) are small compared to the dimension of @xmath31 . the application of gaussian elimination to the extended sparse system of leads to the creation of dense fill - in blocks , which jeopardizes the computational efficiency of the method . in order to obtain a fast solver , it is crucial to preserve the sparsity pattern throughout the elimination . this can be achieved by compressing and redirecting fill - ins that correspond to well - separated interactions , as these are expected to have a numerically low - rank . for example , consider the fill - ins @xmath72 and @xmath73 arising between the well - separated variables @xmath74 and @xmath75 ( the prime indicates a fill - in ; @xmath76 and @xmath77 are originally zero blocks ) . these fill - ins can be approximated by a low - rank representation ( compression step ) by means of a truncated singular value decomposition ( svd ) in which only the most significant singular values and vectors are retained : @xmath78 the rank of approximations and is determined by a prescribed tolerance  @xmath1 ( either relative or absolute ) on the singular values . next , a recompression step is performed to obtain new low - rank interpolation and anterpolation operators @xmath79 and @xmath80 , respectively , as well as @xmath81 and @xmath82 , such that the fill - ins @xmath72 and @xmath73 can be redirected through the existing low - rank interaction between @xmath74 and @xmath75 in @xmath65 . as a result , there is no need to store @xmath72 and @xmath73 explicitly , implying that the sparsity pattern of the matrix is maintained throughout the elimination . if we make the assumption that the rank remains bounded as the problem size increases , it can be demonstrated that this leads to an algorithm with an asymptotic complexity of @xmath9 . although there is no mathematical proof to support this assumption , numerical examples confirm the fact that this is indeed the case in many practical examples  @xcite . the compression and recompression procedures make the ifmm inexact , but the accuracy can be tuned by varying @xmath1 . in terms of computational efficiency , one needs to make a trade - off between a highly accurate direct solver or a low - accuracy preconditioner . we focus in this paper on its use as a preconditioner . more details on the aforementioned procedures are provided in  @xcite . the value of @xmath1 should be chosen in relation to the accuracy of the initial low - rank operators in ( i.e. , the higher the initial rank , the smaller @xmath1 should be ) . several numerical benchmarks are considered in this section to demonstrate the efficiency of the ifmm as a preconditioner for the problems under concern . in all of these benchmarks , is solved for the density vector @xmath83 with a fast multipole accelerated gmres solver  @xcite , and is subsequently used to evaluate the velocities @xmath84 inside the domain . a tolerance of @xmath85 is specified for the relative residual in gmres . the computations have been performed on intel^^ xeon^^ e5 - 2650 v2 ( 2.60  ghz ) cpus .    in all benchmarks , a black - box approach is used in the ifmm for constructing the initial low - rank operators in . more precisely , interpolation based on chebyshev polynomials @xcite is employed ( with @xmath86 chebyshev nodes in each direction ) , followed by an additional svd to reduce the rank . a uniform quadtree decomposition of the domain is used , and the number of levels is adjusted to provide a reasonable trade - off between the time spent at the leaf level and the time spent at higher levels in the tree . we are interested in the porous geometry illustrated in . the outer wall is a rounded off version of the rectangle @xmath87 \times [ -2.6,2.6]$ ] , and the radii of the pores are distributed in @xmath88 $ ] . the smallest gap between two pores is @xmath89 , and the smallest gap between a pore and the outer wall is @xmath90 . finally , the porous region has a porosity of @xmath91 . we test the preconditioner on simpler geometries by considering only the first @xmath92 and @xmath93 left - most pores . we also scale the outer boundary so that the area to the right of the pores remains roughly constant . before studying the pipe flow boundary conditions defined in , we first validate our method by examining the convergence with two sets of boundary conditions that have a closed - form solution in and  [ subsec : analytical ] . we subsequently discuss the pipe flow boundary conditions of in . we use a shear flow boundary condition on all the boundaries so that the exact solution is @xmath94 . with this boundary condition , the density function @xmath83 is only non - zero along the outer boundary @xmath10 . a first small example involves only the first 22 pores ( see ) ; the geometry of the outer boundary @xmath96 is scaled so that the area of the pore - free region is the same as the area of the pore - free region in . we discretize each pore with @xmath97 points and the outer boundary with @xmath98 points . the resulting dense system has nearly @xmath99 unknowns . three different preconditioning approaches are used in gmres : ( i )  no preconditioner , ( ii )  a block - diagonal preconditioner ( with the diagonal blocks corresponding to the self - interactions of the individual pores and the outer boundary ) , and ( iii )  the ifmm . in the ifmm , ten chebyshev nodes are employed in each direction to construct the initial low - rank approximations ( @xmath100 ) , along with a relative accuracy @xmath101 for the low - rank ( re)compressions . the preconditioners in cases ( ii ) and ( iii ) are applied from the left as @xmath102 where @xmath103 represents the preconditioner . \(a ) shows the relative residual as a function of the number of gmres iterations and ( b ) shows the cpu time required by the preconditioner and gmres for all three preconditioning strategies . convergence is very slow without a preconditioner and the relative residual does not drop below the desired tolerance within the prescribed maximum number of 1000 iterations . using the block - diagonal preconditioner , on the other hand , leads to convergence within 139 iterations . this is a very cheap preconditioning strategy , as the time needed to construct the preconditioner is negligible . the ifmm preconditioner is more expensive to compute , but the overall wall time is significantly reduced because the total number of iterations is reduced to merely 40 . note that the relative residual shown in ( a ) is the left preconditioned residual @xmath104 , where @xmath105 is the estimated solution at each gmres iteration . the actual residual @xmath106 is computed as well once the latter drops below the prescribed tolerance . indicates that the discrepancy between the preconditioned and the actual relative residual is very small in all cases . the overall computation times for each of the preconditioning approaches are presented in ( b ) and , where a decomposition is made into the time required to construct the preconditioner and the actual iteration time . it is clear that the additional cost for constructing the ifmm pays off since the total time is more than halved with respect to the block - diagonal preconditioner . ( a ) pores and @xmath107 unknowns . the number of gmres iterations is restricted to 1000.,title="fig : " ] ( b ) pores and @xmath107 unknowns . the number of gmres iterations is restricted to 1000.,title="fig : " ]    .overview of the results for various preconditioning strategies for a shear flow problem with @xmath95 pores and @xmath107 unknowns . the number of gmres iterations is restricted to 1000 . [ cols="^,^,^,^,^ " , ]     shows a quiver plot of the velocity field @xmath18 obtained with the ifmm preconditioner . a very complex flow is resolved in this case as well . obtained when using the ifmm preconditioner for a pipe flow problem with @xmath108 pores and @xmath109 unknowns . the velocity field is plotted in a small window in the middle of the porous region . ] discretization of a first - kind boundary integral equation , representing a 2d incompressible stokes flow in a porous medium , leads to a dense linear system of equations that is computationally expensive to solve . in particular , a large number of iterations is often required if an iterative krylov methods such as gmres is employed . in this paper , the ifmm has been presented as an efficient preconditioner that is capable of significantly reducing the number of iterations and the overall computation cost . the ifmm is in essence an inexact fast direct solver for dense @xmath0-matrices with a linear complexity . this complexity is achieved through two key ideas : the hierarchical low - rank structure of the matrix is exploited to represent the original system as an extended sparse system , and low - rank ( re)compressions are applied to maintain the sparsity pattern of the latter throughout the gaussian elimination . the solver has a tunable accuracy @xmath1 allowing for a trade - off between a highly accurate direct solver and a low - accurate preconditioner . various numerical benchmarks have been carried out to validate the accuracy and to assess the efficacy of the ifmm as a preconditioner . it has been demonstrated that the ifmm preconditioner outperforms a block - diagonal preconditioner for several types of boundary conditions . this is especially the case for pipe flow problems involving many pores , which is relevant for various practical applications .    finally , it is important to note that the current ifmm implementation is based on a uniform quadtree decomposition of the domain , which is far from optimal for the complex geometries considered in this paper . the use of an adaptive quadtree is expected to make the ifmm even more efficient . pc is a post - doctoral fellow of the research foundation flanders ( fwo ) and a francqui foundation fellow of the belgian american educational foundation ( baef ) . the financial support is gratefully acknowledged . bq acknowledges startup funds from florida state university . kitamura k , jiang f , valocchi aj , chiyonobu s , tsuji t , christensen kt . the study of heterogeneous two - phase flow around small - scale heterogeneity in porous sandstone by measured elastic wave velocities and lattice boltzmann method simulation . _ journal of geophysical research : solid earth _ 2014 ; * 119*(10):75647577 . wu d , liu h , xie m , liu h , sun w. experimental investigation on low velocity filtration combustion in porous packed bed using gaseous and liquid fuels . _ experimental thermal and fluid science _ 2012 ; * 36*:169177 . kazemifar f , blois g , kyritsis dc , christensen kt . quantifying the flow dynamics of supercritical co2water displacement in a 2d porous micromodel using fluorescent microscopy and microscopic piv . _ advances in water resources _ 2015 ; . lebon l , oger l , leblond j , hulin jp , martys ns , schwartz lm . pulsed gradient nmr measurements and numerical simulation of flow velocity distribution in sphere packings . _ physics of fluids _ 1996 ; * 8*(2):293301 . de  anna p , borgne tl , dentz m , tartakovsky am , bolster d , davy p. flow intermittency , dispersion , and correlated continuous time random walks in porous media . _ physical review letters _ 2013 ; * 110*(18):184502 . malholtra d , gholami a , biros g. a volume integral equation stokes solver for problems with variable coefficients . _ proceedings of the international conference for high performance computing , networking , storage and analysis _ , ieee press , 2014 ; 92102 . heuer n , stephan ep , tran t. multilevel additive schwarz method for the h - p version of the galerkin boundary element method . _ mathematics of computation of the american mathematical society _ 1998 ; * 67*(222):501518 . nabors k , korsmeyer ft , leighton ft , white jk . preconditioned , adaptive , multipole - accelerated iterative methods for three - dimensional first - kind integral equations of potential theory . _ siam journal on scientific and statistical computing _ 1994 ; * 15*:713735 . langer u , pusch d , reitzinger s. efficient preconditioners for boundary element matrices based on grey - box algebraic multigrid methods . _ international journal for numerical methods in engineering _ 2003 ; * 58*(13):19371953 . sheng z , dewilde p , chandrasekaran s. algorithms to solve hierarchically semi - separable systems . _ system theory , the schur algorithm and multidimensional analysis _ , _ operator theory : advances and applications _ , 176 , alpay d , vinnikov v ( eds . ) . birkhuser basel , 2007 ; 255294 . gillman a , martinsson pg . an @xmath113 algorithm for constructing the solution operator to 2d elliptic boundary value problems in the absence of body loads . _ advances in computational mathematics _ 2014 ; * 40*(4):773796 . coulier p , darve e. efficient mesh deformation based on radial basis function interpolation by means of the inverse fast multipole method . _ computer methods in applied mechanics and engineering _ 2016 ; .
we consider an efficient preconditioner for boundary integral equation ( bie ) formulations of the two - dimensional stokes equations in porous media . while bies are well - suited for resolving the complex porous geometry , they lead to a dense linear system of equations that is computationally expensive to solve for large problems . 
 this expense is further amplified when a significant number of iterations is required in an iterative krylov solver such as gmres . in this paper 
 , we apply a fast inexact direct solver , the inverse fast multipole method ( ifmm ) , as an efficient preconditioner for gmres . 
 this solver is based on the framework of @xmath0-matrices and uses low - rank compressions to approximate certain matrix blocks . 
 it has a tunable accuracy @xmath1 and a computational cost that scales as @xmath2 . 
 we discuss various numerical benchmarks that validate the accuracy and confirm the efficiency of the proposed method . 
 we demonstrate with several types of boundary conditions that the preconditioner is capable of significantly accelerating the convergence of gmres when compared to a simple block - diagonal preconditioner , especially for pipe flow problems involving many pores .