A Hybrid Semi-Lagrangian Cut Cell Method for Advection-Diffusion Problems with Robin Boundary Conditions in Moving Domains
aa r X i v : . [ m a t h . NA ] F e b A Hybrid Semi-Lagrangian Cut Cell Method forAdvection-Diffusion Problems with Robin Boundary Conditions inMoving Domains
Aaron Barrett , Aaron L. Fogelson , and Boyce E. Griffith Departments of Mathematics and Bioengineering, University of Utah, Salt Lake City, UT, USA Departments of Mathematics, Applied Physical Sciences, and Biomedical Engineering, University of North Carolina, ChapelHill, NC, USA Carolina Center for Interdisciplinary Applied Mathematics, University of North Carolina, Chapel Hill, NC, USA McAllister Heart Institute, University of North Carolina, Chapel Hill, NC, USA * [email protected] Abstract
We present a new discretization for advection-diffusion problems with Robin boundary conditions oncomplex time-dependent domains. The method is based on second order cut cell finite volume methodsintroduced by Bochkov et al. [8] to discretize the Laplace operator and Robin boundary condition. Toovercome the small cell problem, we use a splitting scheme that uses a semi-Lagrangian method to treatadvection. We demonstrate second order accuracy in the L , L , and L ∞ norms for both analytic testproblems and numerical convergence studies. We also demonstrate the ability of the scheme to handleconversion of one concentration field to another across a moving boundary. The convective transport and diffusion of chemical species occurs in a broad range of systems. Manyapplications involve chemical concentrations that evolve within complex, time-dependent regions, includingblood flow and clotting in the cardiovascular system [13, 21], particulate and chemical vapor transport in thelungs [14, 35, 33], and drug absorption in the digestive tract [7, 28, 40]. In some cases, critical interactionsoccur between fluid-phase and structure-bound chemicals. These interactions appear in the model as a Robinboundary condition for the fluid-bound chemical. For example, the interaction between circulating proteinsand membrane-bound proteins play a pivotal role in thrombus formation [13]. Modeling these interactionsbecomes even more challenging as one considers the motion of the flow domain itself.The numerical simulation of PDEs in complex domains has garnered significant attention. Embeddedboundary methods are popular approaches in which a fixed rectangular Cartesian grid is overlayed on thecomplex structure. Embedded boundary methods typically involve altering the PDE to include an additionalsource term that is non-zero only near the boundary. For instance, the immersed boundary method [15, 30]uses an integral transform with a regularized delta function kernel to enforce boundary conditions alongirregular interfaces immersed in a background Cartesian grid. These methods typically are designed withDirichlet boundary conditions, and have to be modified with special interpolation procedures to allow for1ther boundary conditions [41, 27, 9]. Volume penalization [18] methods and diffuse domain [39, 22, 36]methods both introduce phase field models to keep track of the interface. A smoothed version of the interfaceis then used to modify the original PDE to account for the boundary conditions, ultimately yielding anapproach similar to the immersed boundary method. These methods all have the effect of smoothing theboundary of the complex domain over several grid cells. In applications in which boundary interactions areimportant, a regularized surface can lower solution accuracy near the boundary. An alternative approach isthe immersed interface method [20, 23, 37], which derives jump conditions across the interface and then buildsthese jump conditions into the discretized equations. Using these methods to impose boundary conditionsrequires relating the boundary conditions to jump conditions. Many applications have succesfully used theimmersed interface method, but to our knowledge, the jump conditions for arbitrary spatially dependentRobin boundary conditions have not yet been determined. A closely related approach is the ghost fluidmethod [10], which uses an interpolation procedure to build stencils near the interface that incorporate theboundary conditions. This approach can lead to non-symmetric and, in some cases, ill-conditioned systemsthat require specialized solvers [38]. An alternative to these methods is the immersed boundary smoothextension method (IBSE) [34], which can solve PDEs on complex geometries by embedding the geometryin a simpler region and solving the PDE on the extended domain, which now includes both “physical” and“non-physical” subdomains. In the IBSE approach, a body force is incorporated in the “non-physical” regimeto extend the physical solution smoothly outside the domain. This allows for high order accuracy to theboundary, but at the cost of solving an additional multiharmonic problem, with the order depending on thenumber of interface conditions.The approach that we use here to impose boundary conditions along a complex boundary while retaininga Cartesian grid discretization framework is based on a cut-cell finite volume framework. In these flux-basedmethod, fluxes are carefully calculated to account for the portion of cells that are inside the physical domain.This approach allows for accurate solutions along the boundary, but comes at the expense of computing cellgeometries at the interface. This expense can be alleviated through the use of a level set function to trackthe surface interface [25]. Recent work by Helgadottir et al. [16] has demonstrated the ability to imposeDirichlet, Neumann, and Robin boundary conditions for Poisson problems. Cut-cell methods can suffer fromthe so-called “small cell problem,” however, in which cells with small volumes necessitate the use of extremelysmall time step sizes for conditionally stable time stepping schemes. This becomes a serious problem forhyperbolic problems for which accurate, efficient, and unconditionally stable implicit time-stepping schemesare difficult to create. Approaches to alleviate the small cell problem for advective PDEs involve using animplicit method only for the cut-cells [24], merging small cells with their neighbors to effectively create alarger cell [32], or through partitioning fluxes into “shielded” and “unshielded” zones based on the geometryof the cut-cell [19]. These methods have been succesfully utilized in two spatial dimensions, but extendingthem to three spatial dimensions remains a challenge. Recently, cut-cell methods that do not suffer from thesmall cell problem have been introduced to handle moving boundaries [32]; however, a formulation involvingRobin boundary conditions has yet to be developed. The major contribution of this study is the constructionof such a method.Herein, we develop a cut-cell method for advection-diffusion equations on moving domains that allowsfor the imposition of general Robin boundary conditions. To avoid the small cell problem, we introduce asplit cut-cell semi-Lagrangian scheme, in which the diffusive operator is handled using well established finitevolume methods, and the advective operator is treated using a semi-Lagrangian method. The benefit ofsplitting the two operators is two fold. First, in the diffusive solve, the flux from the boundary does notneed to account for the change in location of the boundary. This allows us to leverage recent work by Papacet al. [29] and Bochkov et al. [8] to solve a Poisson-like problem with stationary boundaries. Second, the2emi-Lagrangian scheme for advection has no CFL stability constraint, so the small cell problem is no longeran issue. We demonstrate that this method can accurately resolve concentrations near the boundary forboth Robin and Neumann boundary conditions, including inhomogeneous Robin boundary conditions. Inaddition, we show that this method is able to handle conversion of one concentration into another across theboundary, effectively modelling surface reactions.
We consider the advection and diffusion of a chemical species with concentration q ( x , t ) in an arbitrarydomain Ω t , embedded in a larger rectangular region B . Both q ( x , t ) and Ω t are transported by the samevelocity field u ( x , t ). We assume that q ( x , t ) diffuses with diffusion coefficient D . We consider the particularcase in which q ( x , t ) satisfies Robin boundary conditions on the boundary Γ t of Ω t , so that ∂q ( x , t ) ∂t + ∇ · (cid:16) u ( x , t ) q ( x , t ) − D ∇ q ( x , t ) (cid:17) = f ( x , t ) , x ∈ Ω t , (1a) D ∇ q ( x , t ) · n ( x , t ) + a ( x , t ) q ( x , t ) = g ( x , t ) , x ∈ Γ t , (1b)where n ( x , t ) is the outward unit normal of Γ t and f ( x , t ) is a given forcing function. In the implementation, q ( x , t ) is defined only inside Ω t . The bounding region B is used only to define the cells contained within Ω t .We describe the boundary Γ t using a signed distance function φ ( x , t ) such thatΓ t = { x ∈ B | φ ( x , t ) = 0 } . (2)The signed distance function is passively advected by the prescribed velocity, ∂φ ( x , t ) ∂t + ∇ · ( u ( x , t ) φ ( x , t )) = 0 , x ∈ B . (3)Because φ ( x , t ) will not in general remain a signed distance function under advection by (3), a reinitializationprocedure is used to maintain the signed distance property. This can be achieved by computing the steady-state solution to the Hamilton-Jacobi equation ∂ ˆ φ ( x , τ ) ∂τ + sgn (cid:16) ˆ φ ( x , τ ) (cid:17)(cid:16)(cid:13)(cid:13)(cid:13) ∇ ˆ φ ( x , τ ) (cid:13)(cid:13)(cid:13) − (cid:17) = 0 , (4)ˆ φ ( x ,
0) = φ ( x , t ) , after which φ ( x , t ) is set to the steady state solution. To simplify notation, we describe the method in two spatial dimensions. Extensions to a third spatialdimension are straightforward, with the exception of calculating cut-cell geometries. We reference appropriateliterature for modified calculations needed in three spatial dimensions.We use a splitting scheme to split equation (1a) into a diffusion step, ∂q ( x , t ) ∂t − ∇ · ( D ∇ q ( x , t )) = f ( x , t ) , x ∈ Ω t (5a) D ∇ q ( x , t ) · n ( x , t ) + a ( x , t ) q ( x , t ) = g ( x , t ) , x ∈ Γ t , (5b)3here the domain Ω t remains fixed, followed by an advection step, ∂φ ( x , t ) ∂t + ∇ · ( u φ ( x , t )) = 0 for x ∈ B (6a) ∂q ( x , t ) ∂t + ∇ · ( u q ( x , t )) = 0 for x ∈ Ω t , (6b)in which we evolve both the level set and the function q ( x , t ). We note that because Ω t and q ( x , t ) evolve withthe same velocity, no boundary condition is needed with equation (6). The use of this splitting procedureinccurs an additional error cost. This cost can be reduced to second order temporal accuracy using Strangsplitting which will be described later.We overlay a Cartesian grid on top of B such that B consists of rectangular grid cells c i,j and B = ∪ c i,j with a grid cell spacing of ∆ x = ∆ y = h . The concentration field is approximated at the cell centroid ofeach full or partial cell contained within Ω t .In the following sections, we describe the discretization of equations (5) and (6). In what follows, unlessotherwise noted, q i,j refers to the concentration at the cell centroid of the cell c i,j ∩ Ω t . To solve the diffusion step from equations (5), we employ a cut-cell finite volume method based on theapproaches of Papac et al. [29] and Arias et al. [1]. Integrating equation (5a) over a cell c i,j that is entirelyor partially interior to Ω t and dividing by the volume of the cell, we get1 | c i,j ∩ Ω t | Z c i,j ∩ Ω t ∂q ( x , t ) ∂t d x = 1 | c i,j ∩ Ω t | Z c i,j ∩ Ω t D ∆ q ( x , t )d x . (7)We define Q i,j as the cell average of q ( x , t ) in the cell c i,j ∩ Ω t . Replacing the cell average in the left side ofequation (7) and employing the divergence theorem on the right-hand side, we get dQ i,j dt = 1 | c i,j ∩ Ω t | Z ∂ ( c i,j ∩ Ω t ) D ∇ q ( x , t ) n · d A , (8)in which n is the outward unit normal of ∂ ( c i,j ∩ Ω t ). We can further divide the integral in equation (8)into an integral over the cell boundary ∂ c i,j ∩ Ω t and an integral over the physical boundary Γ t ∩ c i,j Z ∂ ( c i,j ∩ Ω t ) D ∇ q ( x , t ) · n dA = Z ∂ c i,j ∩ Ω t + Z c i,j ∩ Γ ! D ∇ q ( x , t ) n · d A . (9)We approximate the first integral by Z ∂ c i,j ∩ Ω t D ∇ q ( x , t ) · n dA ≈ L i + ,j ˆ q i +1 ,j − ˆ q i,j ∆ x − L i − ,j ˆ q i,j − ˆ q i − ,j ∆ x (10)+ L i,j + ˆ q i,j +1 − ˆ q i,j ∆ y − L i,j − ˆ q i,j − ˆ q i,j − ∆ y , in which ˆ q i,j is the point-wise concentration at the center of the cell c i,j and L i + ,j is the length fractionof the face (cid:0) i + (cid:1) × (cid:2) j − , j + (cid:3) covered by the irregular domain (see Figure 1a). It is challenging tocompute L i + ,j exactly. To achieve second order accuracy, it suffices to use a linear approximation L i + ,j = ∆ y (cid:12)(cid:12)(cid:12)(cid:12) φ i + 12 ,j − φ i + 12 ,j − − φ i + 12 ,j + 12 (cid:12)(cid:12)(cid:12)(cid:12) if φ i + ,j − · φ i + ,j + < y if φ i + ,j − < φ i + ,j + <
00 otherwise. (11)4Sfrag replacements ΓΩ t c i,j L i − ,j L i,j − ˆ q i,j +1 ˆ q i,j ∆ yL Γ r i,j d i,j (a) PSfrag replacements ΓΩ t c i,j L i − ,j L i,j − ˆ q i,j +1 ˆ q i,j ∆ y L Γ r i,j d i,j (b) Figure 1: Depiction of the nomenclature used in the diffusion discretization for the cell length fractions (a)and for the boundary fluxes (b).While in the evolution of the level set φ ( x , t ), the degrees of freedom live at cell centers, in equation (11),we require the value of φ ( x , t ) at cell nodes. In our computations, we use a simple bi-linear interpolant tofind these nodal values. In three spatial dimensions, evaluating the corresponding quantity L i + ,j,k wouldinvolve computing the surface area.We remark that ˆ q i,j refers to the cell center of the (possibly cut) grid cell c i,j regardless of the locationof the physical boundary Γ t . In cases where the cell center does not correspond to the location of degreesof freedom, e.g. near cut-cells, we must reconstruct these values. Here, we perform this reconstruction usingeither a moving least squares (MLS) approximation, or a radial basis function (RBF) interpolant. Bothprocedures are described in more detail in Section 3.3.As done in [8], the second integral is approximated using a linear approximation to q ( x , t ) on the boundaryin the direction normal to Γ. Z c i,j ∩ Γ t D ∇ q ( x , t ) · n dA = Z c i,j ∩ Γ t ( g ( x , t ) − aq ( x , t )) dA ≈ L Γ ( g ( r i,j , t ) − aq ( r i,j , t )) , (12)in which r i,j is the closest point on Γ t to x i,j and L Γ = | c i,j ∩ Γ t | (see Figure 1b. The value q ( r i,j ) is foundusing a Taylor series expansion q ( x i,j , t ) = q ( r i,j , t ) + d i,j ∂q ( r i,j , t ) ∂ n + O (cid:0) h (cid:1) , (13)in which d i,j = φ ( x i,j ,t ) |∇ φ ( x i,j ,t ) | is the distance between r i,j and the cell center x i,j . On the boundary, we havethat D ∂q ( r i,j , t ) ∂ n + aq ( r i,j , t ) = g ( r i,j , t ) . (14)We solve equations (13), dropping the O (cid:0) h (cid:1) terms, and (14) for q ( r i,j , t ) and use this value in the dis-cretization of the integral in equation (12). As shown in [8], this system is well posed for a sufficiently refinedgrid.The cut-cell volume | c i,j ∩ Ω t | and physical boundary | c i,j ∩ Γ t | are found by decomposing the cut-celland boundary into simplices, for which analytic formulas for the volume exist [25].We discretize equation (5) using the implicit trapezoidal rule:5∆ t Q n +1 i,j − D | c i,j ∩ Ω − | L i + ,j q n +1 i +1 ,j − q n +1 i,j ∆ x − L i − ,j q n +1 i,j − q n +1 i − ,j ∆ x + L i,j + q n +1 i,j +1 − q n +1 i,j ∆ y − L i,j − q n +1 i,j − q n +1 i,j − ∆ y ! = 1∆ t Q ni,j + D | c i,j ∩ Ω − | (cid:18) L i + ,j q ni +1 ,j − q ni,j ∆ x − L i − ,j q ni,j − q ni − ,j ∆ x (15)+ L i,j + q ni,j +1 − q ni,j ∆ y − L i,j − q ni,j − q ni,j − ∆ y (cid:19) + 1 | c i,j ∩ Ω − | f n + i,j . In our computational examples, we use GMRES without preconditioning to solve this system of equations.The solver typically converges to a solution after approximatley 20 to 50 iterations.
For the advection step, we use a semi-Lagrangian method [11, 31] to advance both the level set φ ( x , t ) andthe concentration q ( x , t ). We solve equation (6) by first finding the preimage X i,j = χ − (cid:0) x i,j , t n +1 (cid:1) of thefluid parcel located at x i,j at time t n +1 . The mapping χ ( x , t ) satisfies the differential equation ∂ χ ( x , t ) ∂t = u ( x , t ) for x ∈ B . (16)The preimage X i,j can be found by integrating equation (16) backwards in time. In this work, we use anexplicit two step Runge Kutta method: X ⋆i,j = x i,j − ∆ t u (cid:0) x i,j , t n +1 (cid:1) , (17a) X ni,j = x i,j − ∆ t u (cid:16) X ⋆i,j , t n + (cid:17) . (17b)With the preimage found, we can interpolate the solution at time n to the preimage locations X ni,j .The interpolation procedure consists of one of two methods depending on whether X ni,j is near cut-cells.Away from cells pierced by the zero level set, we use a tensor product of special piecewise Hermite polynomialscalled Z-splines [6]. The Z-spline Z m ( x ) interpolating the data ( x i , f i = f ( x i )) ni =1 is a piecewise polynomialfunction that satisfies Z m ( x ) ∈ C m ([ x , x n ]) , (18a) d p dx p Z m ( x ) (cid:12)(cid:12)(cid:12)(cid:12) x j = f pm,j for p = 0 , . . . , m and j = 1 , . . . , n, (18b) Z m ( x ) ∈ π m +1 ([ x i , x i +1 ]) for i = 1 , . . . , n − , (18c)in which f pm,j is the approximation of the p th order derivative of f ( x ) computed from high-order finitedifferences of f j using 2 m + 1 points, and π m +1 is the space of polynomials of degree less than or equal to2 m + 1.It is possible to define Z m for a general set of data points (or abscissae) via the cardinal Z-splines [6] Z m ( x ) = X i f i ˜ Z m ( x − x i ) . (19)The cardinal Z-splines have two key properties. First, the cardianl Z-splines have compact support, ˜ Z m ( x ) =0 for | x | > m + 1, which allows for efficient local evaluation of interpolants. Second, cardinal Z-splines are6nterpolatory, which makes the interpolant trivial to form. Because data values are defined on a regularCartesian grid, the computation of the full interpolant uses a tensor product of cardinal Z-splines Z m ( x ) = X i X j q i,j ˜ Z m ( x − x i ) ˜ Z m ( y − y j ) . (20)In our computations, we use the quintic Z-splines, which are defined in terms of Z ( x ) = − x − x + x − x if | x | ≤ − x − x + x − x + x if 1 < | x | ≤ − x + x − x + x − x if 2 < | x | ≤
30 otherwise. (21)In cases where we do not have enough points to define the Z-spline, e.g., near cut-cells, we use theprocedure described in Section 3.3. We note that although the interpolation procedure described here has adiscontinuous switch between operators, this does not affect the overall convergence rates.The evolution of the level set and the concentration use the same procedure, except for the choice ofreference grid. The reference grid for the level set consists of the entire domain B , whereas the referencegrid for the concentration consists of the domain Ω t n +1 = (cid:8) x ∈ B | φ (cid:0) x , t n +1 (cid:1) < (cid:9) . Therefore, we updatethe level set prior to the update of the concentration. Near cut-cells, it is necessary to form interpolants on unstructured data. In this study, we use either a movingleast squares (MLS) approximation, or a local radial basis function (RBF) interpolant. This procedure isused in the diffusion step to extrapolate data from cut-cell centroids to full-cell centers and in the advectionstep for cells where the Z-spline interpolant can not be formed. We describe the procedure in the context ofreconstructing a function f ( x ) from the data points { x i , f ( x i ) } with i = 1 , , . . . , n where n is an arbitrarynumber of points.The MLS approximation is formed by finding the best approximation q ( x ) at a point x c of the form q ( x ) = n X j =1 c j p j ( x ) (22)in which the p j ( x ) form a basis for the space of polynomials. The approximation q ( x ) is chosen with respectto the standard weighted L inner product with weight function w ( x ) h f ( x ) , g ( x ) i = Z Ω f ( x ) g ( x ) w ( x ) d x . (23)We use a stencil width of approximately two grid cells. It has been shown that the weight function must besingular at the data locations x i for the MLS approximation to interpolate the data [5]. Here, we use w ( x ) = e −k x − x c k . (24)Because this weighting function is not singular at x c , the reconstructed polynomial will not be interpolatory.The MLS calculation reduces to solving a linear system of the formΛ A c = Λ f (25)in which A is a matrix whose entries consist of a i,j = p j ( x i ) and Λ is a diagonal matrix consisting of elements λ i,i = w ( x i ). In our computational examples, we have observed that using a lower order reconstruction for7he diffusion step is more robust without affecting the order of accuracy. Therefore, we use a quadraticpolynomial in the reconstructions for the advection step and a linear reconstruction with the diffusion step.The RBF approximation is constructed via a polyharmonic spline of the form q ( x ) = n X j =1 λ j φ ( k x − x j k ) + s X j =1 β j p j ( x ) (26)in which φ ( x ) is a polyharmonic radial basis function of degree m and p j ( x ) are a set of s polynomial basisfunctions [5, 12]. The degree of the polynomials k is chosen such that m = 2 k + 1. The coefficients arechosen so that q ( x j ) = f j for j = 1 , . . . , n, (27) n X k =1 λ j p k ( x j ) = 0 for j = 1 , . . . , n. (28)This leads to the linear system for the coefficients A PP T ! λβ ! = f ! (29)in which A is a square matrix with elements A i,j = φ ( k x i − x j k ), P is a rectangular matrix with elements P i,j = p j ( x i ). In our computations, we use the polyharmonic spline φ ( x ) = k x k combined with linearpolynomials. We summarize the full procedure to advance the solution from time t n to time t n +1 :1) If needed, reinitialize the level set φ ( x , t ) into the signed distance function by iterating equation (4) tosteady state using the procedure as described by Nangia et al. [26].2) Update the diffusion equation to half time t n + using the methods described in Section 3.1.3) Update the level set using the prescribed velocity u to time t n +1 by the methods described in Section3.2.4) Advect the concentration using the prescribed velocity u by the methods described in Section 3.2.5) Update the diffusion equation to full time t n +1 using the concentration and level sets from the previoustwo steps.The above procedure is implemented using the SAMRAI [17] infrastructure which provides an efficient,parallelized environment for structured adaptive mesh refinement. The diffusion solve is computed usingsolvers provided by PETSc [4, 2, 3]. Here we demonstrate the capabilities of the method both using a prescribed level set and an level setadvected with the fluid. We start with diffusion dominated examples before exploring inclusion of advectionand spatially varying boundary conditions. 8 .1 Diffusion from a point source
We consider the advection-diffusion of a point source within a disk of radius R . The disk is advected withvelocity u = (cid:0) cos (cid:0) π (cid:1) , sin (cid:0) π (cid:1)(cid:1) T . The exact concentration is q ( x , t ) = 104 D (cid:0) t + (cid:1) e − k x − x c( t ) k D ( t + 12 ) for x ∈ B , (30)in which x c ( t ) is the center of the disk. We apply Robin boundary conditions of the form D ∂q∂ n = g ( x , t ) − a ( x , t ) q ( x , t ) . (31)We set a ( x , t ) = 1 and use the method of manufactured solutions to determine the value of g ( x , t ). Forthis example, we specify the location of the level set φ at each timestep. This procedure allows us test theconvergence of the advection-diffusion step without having to account for numerical error in the evolution ofthe level set. We note that whereas this example includes a trivial semi-Lagrangian backwards integrationstep, the interpolation procedure and diffusion step are non-trivial.The extended domain B consists of the box [0 , × [0 , N points ineach direction. We use a disk radius of R = 1 and initial center x c (0) = (1 . , . T . The choice of centerslightly offsets the disk from the background grid so that different cut-cells are generated on opposite sidesof the disk. We use a diffusion coefficient of D = 0 .
01. The simulations are run at an advective CFL numberof C = ∆ x | u | ∆ t = 0 . L , L , and L ∞ norms.The coarsest simulation uses N = 128 grid points in each direction of the background Cartesian mesh, whichcorresponds to approximately 26 grid cells covering the diameter of the disk. Despite the relatively coarsedescription of the disk, we still see errors that are on the order of one percent.We also consider the case where the level set φ is advected with the same velocity. In this case, numericalissues make it difficult to assess the pointwise error on the Cartesian grid because the piecewise linearreconstruction used to determine the cut-cells gives different cut-cell geometries than that of a prescribedlevel set. Instead, we compare the solution at the boundary to the imposed boundary condition. Wecompute the solution on the boundary using a moving least squares linear extrapolation from cut-cells tothe cell boundary. Further, to assess the convergence rate, we compute the error of the surface integral E = (cid:12)(cid:12)(cid:12)(cid:12)Z π q ( x , t ) dθ − Z π ˜ q ( x , t ) dθ (cid:12)(cid:12)(cid:12)(cid:12) , (32)in which θ is the angle along the boundary from the center of the disk, q ( x , t ) is the exact solution from (30),and ˜ q ( x , t ) is the approximate solution computed using the moving least squares extrapolation as described inSection 3.3. A convergence plot is shown in Figure 2 and indicates that the method achieves between first andsecond order convergence rates. In both the evolved and prescribed level set cases, the RBF reconstructionshows errors that are more than an order of magnitude smaller than the MLS reconstructions. In this section, we again consider diffusion from a point source, but instead apply a solid body rotation tothe disk. In this case, the semi-Lagrangian sub-step contains a non-trivial integration in time, and thereforetests all components of the solver. We again prescribe the initial condition as in equation (30), and applyboundary conditions as in equation (31). 9 .10.060.040.020.01 Δx −4 −3 E rr o r RBFL ΔrateΔ=Δ2.095L ΔrateΔ=Δ2.104L ∞ ΔrateΔ=Δ2.144 (a) Δx −2 −1 E rr o r MLSL ΔrateΔ=Δ1.884L ΔrateΔ=Δ1.886L ∞ ΔrateΔ=Δ1.914 (b) Δx −4 −3 −2 −1 E rr o r ErrorΔonΔBoundaryRBFΔRateΔ=Δ1.745MLSΔRateΔ=Δ1.812 (c)
Figure 2: Convergence rates for advection diffusion from a point source in a constant flow using either aRBF (a) or MLS (b) reconstruction. Also shown is the convergence rate for the integral in equation (32)using both an RBF or MLS reconstruction (c). In all cases, RBF reconstructions are greater than an orderof magnitude more accurate than MLS reconstructions.The computational domain B is [ − , × [ − , R = 1 and initial center x c (0) = (1 . , . T . We use a diffusion coefficient of D = 0 .
1. The prescribed velocity is u = 2 π ( − y, x ) T .Simulations are run at an advective CFL number of 0 . T = 1. Figure 3 shows the initialand final solution. We can see that the solution maintains symmetry throughout the simulation. Figure 4shows convergence rates for both RBF and MLS reconstructions. Figures 4a and 4b show convergence ratesin which the level set is prescribed. We see second order convergence rates in all norms for the prescribedlevel set case. Figure 4c shows the convergence rate for the solution extrapolated to the boundary for thecase in which the level set is evolved with the fluid velocity. In this case, we see approximately secondorder convergence rates. In both of these cases, the error for the RBF reconstruction is almost an order ofmagnitude lower than that of the MLS reconstruction. We now consider the advection and diffusion of a concentration inside a vesicle under oscillating Couetteflow with no flux boundary conditions. Specifically, we specify the rotational component of the velocity fieldas u θ ( r ) = (cid:18) ar + br (cid:19) sin( πt ) , (33a) a = Ω R − Ω R R − R , (33b) b = (Ω − Ω ) R R R − R , (33c)where R and R are the radii of the coaxial cylinders, and Ω and Ω are the angular velocities of thecylinders. Here, we set the radii to be R = 0 . R = 3 .
75 with respective angular velocities Ω = 10 . a) (b) Figure 3: Initial (a) and final (b) configuration for diffusion from a point source in a rigid body rotationalflow field. The zero contour of φ is shown in red. Δx −1 E rr o r RBFL ΔrateΔ=Δ2.004L ΔrateΔ=Δ1.998L ∞ ΔrateΔ=Δ2.030 (a) Δx E rr o r MLSL ΔrateΔ=Δ1.940L ΔrateΔ=Δ1.878L ∞ ΔrateΔ=Δ1.533 (b) Δx E rr o r ErrorΔonΔBoundaryRBFΔRateΔ=Δ1.944MLSΔRateΔ=Δ1.490 (c)
Figure 4: Convergence rates for advection diffusion from a point source in a solid body rotational flow usingeither a RBF (a) or MLS (b) reconstruction. Also shown is the convergence rate for the integral in equation(32) using both an RBF or MLS reconstruction (c). 11nd Ω = 1 . φ v ( x , t ) so that the initial condition is φ v ( x ,
0) = k x − x v , c k − R (34)where R = 1 is the radius of the disk, and x v , c = (1 . , . T is the center of the disk. We initialize theconcentration field q v via q v ( x ,
0) = (cos( π k x − x v , c k ) + 1) (35)For comparison, we define another concentration field q o ( x , t ) outside the vesicle, but inside the cylinders.The second level set is defined by φ o ( x , t ) = max( − φ v ( x , t ) , k x k − R , R − k x k ) . (36)The initial condition is given by q o ( x ,
0) = ( (cos( π k x − x o , c k ) + 1) if k x − x o , c k ≤
10 otherwise, (37)in which x o , c = − x v , c .The computational domain consists of the box [ − , × [ − , D = 0 .
05. Figure 5 shows the solution at five different timevalues during the simulation. Over the course of the simulations, the concentration field inside the disk stayscompletely contained within the disk, while the outside concentration field is free to diffuse between thecoaxial cylinders.Figures 6 and 7 show numerical convergence results for the concentrations inside the vesicle and outsidethe vesicle but inside the two cylinders respectively. We estimate the convergence rate r using Richardsonextrapolation via r = log k q h − q h kk q h − q h k , (38)in which q h is the solution with a grid spacing of h . To compute the difference k q h − q h k , we first interpolatethe fine solution q h onto the coarser grid q h . Because the interpolation from a fine grid to a coarse grid isnontrivial near cut-cells, we exclude all cut-cells and any cells that neighbor cut-cells in the computation ofthe norm. To assess the error on the cut-cells, we perform a reconstruction of the the solution on the surfaceof the disk by extrapolating the solution to the disk in each cut-cell, then computing a cubic spline of thesedata points. Figure 6d shows convergence of this reconstruction for the solution inside the vessicle. Weobtain second order convergence rates using RBFs at all points in time for the finest grids. In contrast, theMLS reconstruction does not show optimal convergence rates for the range of grid spacings considered. Weexpect that under additional grid refinement the MLS approximation will settle to second order accuracy,but this is not achieved over the range of grid spacings considered herein. The dips in convergence rates inFigure 7 occur when significant values of q o begin to appear in cut cells. Additionally, Figures 8 and 9 showthe norms of the differences. To more fully qualify the difference between MLS and RBF reconstructions,we estimate the error coefficient C in the expression k q h − q k ≈ Ch p , (39)in which q is the exact solution and p is the order of accuracy. We approximate C by comparing differencesbetween two grids of refinement h and 2 h through the following equation C ≈ k q h − q h k h p (1 − p ) . (40)12 o q v L norm L norm L ∞ norm L norm L norm L ∞ norm C MLS C RBF C coefficient for both MLS and RBF reconstructions. In all cases, the RBF reconstructionyields improved accuracy.For both MLS and RBF reconstructions, we have p = 2. Table 1 shows our estimates of the C coefficient forboth reconstructions. In all cases, the RBF reconstruction is more accurate. We now consider the case of two fluid phase chemicals advecting and diffusing on two different domains, butinteracting through a common boundary. Specifically, we consider oscillatory Couette flow of two chemicals,one of which is contained within a vesicle that passively advects with the flow and is initially in the shape ofa disk. The other chemical exists outside the vesicle, but between the disks defining the domain of interest.Specifically, we solve the equations ∂q v ∂t + u · ∇ q v = D ∆ q v , x ∈ Ω v (41a) − D ∂q v ∂ n = κ ( q v − q o ) , x ∈ Γ v (41b) ∂q o ∂t + u · ∇ q o = D ∆ q o , x ∈ Ω o (41c) − D ∂q o ∂ n = κ ( q o − q v ) , x ∈ Γ v (41d) − D ∂q o ∂ n = 0 , x ∈ Γ o \ Γ v , (41e)in which u is defined in the previous section and Ω v and Ω o are the domains for the concentration insideand outside the vesicle, respectively. The initial concentration for q v is given by equation (35), and for q o isinitially uniformly zero.We modify the time discretization of the diffusion step to use a modified trapezoidal rule for the boundaryconditions. We first solve equation (15) for both q n +1o and q n +1v using explicit approximations for theboundary conditions − D ∂q n +1o ∂ n = κ (cid:0) q n +1o − q n v (cid:1) (42) − D ∂q n +1v ∂ n = κ (cid:0) q n +1v − q n o (cid:1) . (43)We then solve equation (15) again using the intermediate result when evaluating the boundary conditions − D ∂q n +1o ∂ n = κ (cid:0) q n +1o − ˜ q v (cid:1) (44) − D ∂q n +1v ∂ n = κ (cid:0) q n +1v − ˜ q o (cid:1) , (45)13Sfrag replacements4 . . . . . . . . . (a) PSfrag replacements4 . . . . . . . . . (b) PSfrag replacements4 . . . . . . . . . (c) PSfrag replacements4 . . . . . . . . . (d) PSfrag replacements4 . . . . . . . . . (e) PSfrag replacements 4 . . . . . . . . .
30 PSfrag replacements4 . . . .
55 0 . . . . . q o (orange) and q v (blue) at time points 0 . . . . . φ o is shown in green.14 .00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00Time1.001.251.501.752.002.252.50 R a t e s Inside L ratesMLS 128-512MLS 256-1024RBF 128-512RBF 256-1024 (a) R a t e s Inside L rates (b) R a t e s Inside L ∞ rates (c) R a t e s Inside Boundary L rates (d) Figure 6: Numerical convergence study for the interior of the disk in the L (a), L (b), and L ∞ (c)norms. The convergence rates are computed from simulations with N = 128 , ,
512 points (blue) and N = 256 , , R a t e s Outside L ratesMLS 128-512MLS 256-1024RBF 128-512RBF 256-1024 (a) R a t e s Outside L rates (b) R a t e s Outside L ∞ rates (c) Figure 7: Numerical convergence study for the exterior of the disk but inside the two cylinders in the L (a), L (b), and L ∞ (c) norms. The convergence study is computed from simulations with N = 128 , , N = 256 , , q o near the boundaries.15 .00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00Time10 −4 −3 −2 D i ff e r e n c e i n n o r m s Inside L normsMLS 128-256MLS 256-512MLS 512-1024RBF 128-256RBF 256-512RBF 512-1024 (a) −4 −3 −2 D i ff e r e n c e i n n o r m s Inside L norms (b) −4 −3 −2 D i ff e r e n c e i n n o r m s Inside L ∞ norms (c) −4 −3 −2 −1 D i ff e r e n c e i n n o r m s Inside Boundary L norms (d) Figure 8: Norms of solution differences for the concentration field inside the domain for the L (a), L (b),and L ∞ (c) norms. The dashed lines use a MLS reconstruction while the solid lines use a RBF reconstruction.Also shown is the norms of solution differences for the extrapolation to the boundary (d). −4 −3 −2 D i ff e r e n c e i n n o r m s O tside L normsMLS 128-256MLS 256-512MLS 512-1024RBF 128-256RBF 256-512RBF 512-1024 (a) −4 −3 −2 D i ff e r e n c e i n n o r m s Outside L norms (b) −4 −3 −2 D i ff e r e n c e i n n o r m s Outside L ∞ norms (c) Figure 9: Norms of solution differences for the concentration field outside the domain for the L (a), L (b),and L ∞ (c) norms. The dashed lines use a MLS reconstruction while the solid lines use a RBF reconstruction.16n which ˜ q v and ˜ q o are the results from the intermediate result.For this example, we perform only RBF reconstructions at the boundaries. We set the parameter κ tobe one. The remaining parameters are the same as in the previous section. Concentration values for variouspoints in time are shown in Figure 10. As q v reaches the boundary Γ v , it gets converted to q o at a rateproportional to the difference in concentrations.As in the previous section, we perform a numerical convergence study for both the interior and exteriorfields. Figures 11 and 12 show the convergence rates as a function of time. The convergence rates as we refinethe grid appear to be converging toward second order, although the method has not yet settled down intoit’s asymptotic regime. We also perform a convergence study for the solution extrapolated to the boundary,which is shown in Figure 11d. We again see approximate second order rates for this reconstruction. A plotof the norms of the solution differences is shown in Figures 13 and 14. The differences in the finest grid showa difference of less than one percent at the final time.To assess conservation, we perform a convergence study on the total amount of the two concentrationfields. Figure 15 shows the total amounts of both fields as a function of time, while the sum of both amountsremains relatively constant. Despite this method not being conservative, we see a loss of less than onepercent over the course of the finest simulation. The total amount of concentration converges at a rate thatis between first and second order accuracy. We have presented a numerical method to simulate advection-diffusion problems with Robin boundaryconditions on irregular, evolving domains. The method shows second order convergence in the L , L ,and L ∞ norms. We have also demonstrated the ability to accurately reconstruct solution values on theboundary, achieving second order accurate results. We have further demonstrated the method to be ableto capture transformation of one concentration into another concentration with all interaction mediatedthrough boundary conditions. Although only two dimensional tests are presented, we expect this method tobe easily extended to three spatial dimensions. The use of radial basis functions when compared to movingleast squares gives optimal convergence rates and more accurate solutions for the grid spacings considered. Z-splines were used in this study for computational efficiency, although any suitable reconstruction procedureshould work. This method has many possible applications, including models with chemical transport inevolving bodies, such as esophageal transport, oxygen flow in the lungs, and blood flow. Acknowledgements
The authors thank Dr. Varun Shankar for helpful discussions on both radial basis functions and semi-Lagrangian schemes. This research was supported through NHLBI award 5U01HL143336 and NSF awardsDMS 1664645, OAC 1450327, and OAC 1931516.
References [1] V. Arias, D. Bochkov, and F. Gibou. “Poisson equations in irregular domains with Robin bound-ary conditions - Solver with second-order accurate gradients”.
Journal of Computational Physics . . . . . . . . . (a) PSfrag replacements0 . . . . . . . . . (b) PSfrag replacements0 . . . . . . . . . (c) PSfrag replacements0 . . . . . . . . . (d) PSfrag replacements0 . . . . . . . . . (e) PSfrag replacements 0 . . . . . . . . . . . . .
075 0 . . . . . q o (orange) and q v (blue) at times 0 . . . . . φ o is shown in green.18 .50 0.75 1.00 1.25 1.50 1.75 2.00Time1.21.41.6 R a t e s Inside L rates 128-512256-1024512-2048 (a) R a t e s Inside L rates (b) R a t e s Inside L ∞ rates (c) R a t e s Inside Boundary L rates (d) Figure 11: Numerical convergence study for the interior of the disk in the L (a), L (b), and L ∞ (c)norms. The convergence rates are computed from simulations with N = 128 , ,
512 points (blue), N =256 , , N = 512 , , R a t e s Outside L rates128-512256-1024512-2048 (a) R a t e s Outside L rates (b) R a t e s Outside L ∞ rates (c) Figure 12: Numerical convergence study for the exterior of the disk but inside the two cylinders in the L (a), L (b), and L ∞ (c) norms. The convergence study is computed from simulations with N = 128 , , N = 256 , , N = 512 , , .0 0.5 1.0 1.5 2.0Time10 −5 −4 −3 −2 −1 D i ff e r e n c e i n n o r m s Inside L norms RBF 128-256RBF 256-512RBF 512-1024RBF 1024-2048 (a) −5 −4 −3 −2 −1 D i ff e r e n c e i n n o r m s Inside L norms (b) −4 −3 −2 −1 D i ff e r e n c e i n n o r m s Inside L ∞ norms (c) −3 −2 −1 D i ff e r e n c e i n n o r m s Inside Boundary L norms (d) Figure 13: Norms of solution differences for the concentration field inside the domain for the L (a), L (b),and L ∞ (c) norms. As shown is the norms of the differences between solutions extrapolated to the boundary(d). −6 −5 −4 −3 −2 −1 D i ff e r e n c e i n n o r m s Outside L norms RBF 128-256RBF 256-512RBF 512-1024RBF 1024-2048 (a) −5 −4 −3 −2 −1 D i ff e r e n c e i n n o r m s Outside L norms (b) −5 −4 −3 −2 −1 D i ff e r e n c e i n n o r m s Outside L ∞ norms (c) Figure 14: Norms of solution differences for the concentration field outside the domain for the L (a), L (b), and L ∞ (c) norms. 20 .0 0.5 1.0 1.5 2.0Time0.00.51.01.52.0 A m o un t Total AmountsInside amountOutside amountTotal amount (a) −2 Δx10 −2 −1 E rr o r ConvergenceΔofΔTotalΔAmountrateΔ=Δ1.614 (b)
Figure 15: Total concentration of both q v and q o (a). The total amount converges at a rate that is betweenfirst and second order (b).[2] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V.Eijkhout, W. Gropp, D. Karpeyev, D. Kaushik, M. Knepley, D. May, L. C. McInnes, R. T. Mills, T.Munson, K. Rupp, P. Sanan, B. Smith, S. Zampini, H. Zhang, and H. Zhang. PETSc Users Manual .Tech. rep. ANL-95/11 - Revision 3.13. Argonne National Laboratory, 2020.[3] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V.Eijkhout, W. Gropp, D. Karpeyev, D. Kaushik, M. Knepley, D. May, L. C. McInnes, R. T. Mills,T. Munson, K. Rupp, P. Sanan, B. Smith, S. Zampini, H. Zhang, and H. Zhang.
PETSc Web page . . 2019.[4] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith. “Efficient Management of Parallelism inObject Oriented Numerical Software Libraries”. Modern Software Tools in Scientific Computing . Ed.by E. Arge, A. M. Bruaset, and H. P. Langtangen. Birkh¨auser Press, 1997, pp. 163–202.[5] V. Bayona. “Comparison of Moving Least Squares and RBF+poly for Interpolation and DerivativeApproximation”.
Journal of Scientific Computing
SIAM Journal on Scientific Computing
Drug Discovery Today
Journal of Computational Physics
376 (2019), pp. 1156–1198.[9] F. Bouchon and G. H. Peichl. “An Immersed Interface Technique for the Numerical Solution of theHeat Equation on a Moving Domain”.
Numerical Mathematics and Advanced Applications 2009 . Berlin,Heidelberg: Springer Berlin Heidelberg, 2010, pp. 181–189.[10] M. Chai, K. Luo, C. Shao, H. Wang, and J. Fan. “A finite difference discretization method for heatand mass transfer with Robin boundary conditions on irregular domains”.
Journal of ComputationalPhysics
400 (2020), p. 108890.[11] D. R. Durran. “Semi-Lagrangian Methods”. 1999, pp. 303–333.[12] N. Flyer, B. Fornberg, V. Bayona, and G. A. Barnett. “On the role of polynomials in RBF-FD ap-proximations: I. Interpolation and accuracy”.
Journal of Computational Physics
321 (2016), pp. 21–38. 2113] A. L. Fogelson and K. B. Neeves. “Fluid Mechanics of Blood Clot Formation”.
Annual Review of FluidMechanics
Toxicological Sciences
Annual Reviewof Fluid Mechanics
Computers and Fluids
121 (2015), pp. 68–80.[17] R. D. Hornung and S. R. Kohn. “Managing application complexity in the SAMRAI object-orientedframework”.
Concurrency and Computation: Practice and Experience
Journal of ComputationalPhysics
Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engi-neering Sciences
Journal of Computational Physics
400 (2020), p. 108854.[21] K. Leiderman and A. L. Fogelson. “Grow with the flow: A spatial-temporal model of platelet depositionand blood coagulation under flow”.
Mathematical Medicine and Biology
Communicationsin Mathematical Sciences
Applied Mathematics and Mechanics
Journal of Scientific Computing
Journal of Computational Physics
Journal of Computational Physics
390 (2019), pp. 548–594.[27] A. Pacheco-Vega, J. R. Pacheco, and T. Rodi´c. “A General Scheme for the Boundary Conditions inConvective and Diffusive Heat Transfer With Immersed Boundary Methods”.
Journal of Heat Transfer
Drug Metabolism and Disposition
Journal of Computational Physics
Acta Numerica
11 (2002), pp. 479–517.[31] G. Russo and F. Filbet. “Semilagrangian schemes applied to moving boundary problems for the BGKmodel of rarefied gas dynamics”.
Kinetic & Related Models
Journal of Computational Physics
235 (2013), pp. 786–809.[33] J. D. Schroeter, J. S. Kimbell, and B. Asgharian. “Analysis of particle deposition in the turbinateand olfactory regions using a human nasal computational fluid dynamics model”.
Journal of AerosolMedicine
Journal of ComputationalPhysics
335 (2017), pp. 155–178.[35] G. Tian and P. W. Longest. “Development of a CFD boundary condition to model transient vaporabsorption in the respiratory airways”.
Journal of Biomechanical Engineering
Journal of Compu-tational Physics
361 (2018), pp. 424–441.[37] S. Xu and Z. J. Wang. “An immersed interface method for simulating the interaction of a fluid withmoving boundaries”.
Journal of Computational Physics
International Journal for NumericalMethods in Fluids
Modelling andSimulation in Materials Science and Engineering
Advanced Drug DeliveryReviews