Space-Time Localized Radial Basis Function Collocation Methods for PDEs Alfa Heryudono Department of Mathematics University of Massachusetts Dartmouth ICERM Topical Workshop Providence August 2017 This research is partly supported by NSF DMS-1552238. Joint work with grad student: Jacob Sousa. Computations are mostly done on UMassD Rapid Prototyping Server. 1 / 21
Dealing with Time-Dependent PDEs for RBF Methods Method of Lines RBF discretization in space + common ODE solver in time. Min changes of PS/FD codes: replace differentiation matrices with RBF versions (Global, RBF-FD, RBF-PU, etc). PS/FD treatments for BCs: Strip-rows, Strip-rows move over columns, fictitious pts/ rect projection (for multiple bcs), penalty, etc. Stability for linear pde case: Eigenvalue and Pseudospectra. Simultaneous Space-Time RBF Boundary value collocation problem in space-time domain. Time is treated as another space variable. RBF-BVP solver have been studied for quite a while. Less worry about choosing ODE solver based on PDE types. Adaptivity, moving boundary, and BCs: same treatments as in BVP cases. No need to rewrite the pde due to var trans (e.g in moving boundary case). Analyzing stability is not clear (e.g. in moving boundary case). Might be expensive to solve (e.g. finding preconditioner, non-linear case). 2 / 21
Space-Time PS Collocation Method: 1D+t linear case PDE : u t = u x ( x , t ) ∈ [ − 1 , 1) × (0 , T ] IC : u ( x , 0) = f ( x ) BC : u (1 , t ) = g ( t ) Use PS or Block PS (Driscoll-Fornberg) to create differentiation matrices. 1 0.75 1 0.5 0.75 1 0.5 0.25 0 0.25 -1 -1 0 -0.5 0 0 -1 -0.5 0 0.5 1 0.5 1 3 / 21
Space-Time PS Collocation Method: 2D+t, linear case PDE : u t = ∆ u + F ( x , y , t ) 2 ( x , y , t ) ∈ Ω × (0 , T ] 1 IC : u ( x , y , 0) = f ( x , y ) BC : u ( ∂ Ω , t ) = g ( ∂ Ω , t ) 0 1 0.5 1 0.5 kron ’s disease is worse in 2 D + t case. 0 0 -0.5 -0.5 -1 -1 10 2 0 P = symrcm(PLinop); L = gpuArray(Linop(P,P)); 1000 10 -2 PL = gpuArray(PLinop(P,P)); 2000 r = gpuArray(rhs(P)); 10 -6 MAXITER = 30; TOL = 1e-14; RESTART = []; 3000 [Ugpu,FLAG,RELRES,ITER,RESVEC] = ... 10 -10 4000 gmres(L,r,RESTART,TOL,MAXITER,PL); U(P) = gather(Ugpu); 0 1000 2000 3000 4000 10 -14 nz = 163249 0 10 20 30 4 / 21
Space-Time PS Collocation Method: 1D+t, nonlinear case Human tear film dynamics: 1D model: see H. et. al 2007 h t + q x = 0 on X ( t ) ≤ x ≤ 1 , where � h 3 � 3 + β h 2 q ( x , t ) = Sh xxx Boundary conditions h ( X ( t ) , t ) = h (1 , t ) = h 0 q ( X ( t ) , t ) = X t h 0 + Q top q (1 , t ) = − Q bot . 0 12 t 8 h(x,t) 4 3.52 0 −1 −0.5 0 0.5 1 x Advance the solution in space-time domain: Slab by Slab (Show MATLAB ). 5 / 21
bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc RBF-FD Differentiation Matrices n loc ∂ Ω x 11 x 27 � λ k φ k ( x ) , s j ( x ) = x 2 x 35 x 8 k =1 ⊕ ⊕ x 33 x 1 x 21 where φ k ( x ) is a radial basis func- Ω x 41 x 18 tion centered at x k . Or in Lagrange formulation as n loc � Ψ k ( x ) u k , s j ( x ) = k =1 where Ψ 1 ( x ) φ 1 ( x ) A − 1 � � Ψ n loc ( x ) � � φ n loc ( x ) � � Ψ = · · · = · · · , � Ψ 1 � φ 1 x ( x ) � � A − 1 � Ψ n loc φ n loc x ( x ) � Ψ x = x ( x ) · · · = x ( x ) · · · , The matrix A with entries a ℓ k = φ k ( x ℓ ) , ℓ, k = 1 , . . . , n loc is local RBF interpolation matrix . BYODM: Bring Your Own Differentiation Matrices 6 / 21
Getting the space-time domain This is probably for programming on a lazy Sunday: Use Mathematica ’s DiscretizeRegion family commands. Surprisingly, Mathematica has many built-in funky domains too. This is also useful if you want to compare results with finite-element. R = ImplicitRegion[-0.6 Sin[t] <= x, {{x, -1, 1}, {t, 0, 1.5 Pi}}]; ev = DiscretizeRegion[R]; pts = MeshCoordinates[ev]; Export["spacetimedom.mat", pts]; To obtained boundary points, you can use Mathematica or boundary command in MATLAB . 7 / 21
t+1D Advection Example solution in space-time domain PDE : u t = u x ( x , t ) ∈ [ X ( t ) , 1) × (0 , T ] IC : u ( x , 0) = f ( x ) BC : u (1 , t ) = g ( t ) 1+( ε r ) 2 . r 2 = ( x − x i ) 2 + ( t − t i ) 2 1 √ IMQ-RBF: f ( x ) = e − 10( x − 0 . 15+0 . 35 y ) 2 , g ( t ) = 0 D t − D x 0 = u f 0 I g P = symrcm(L); u(P) = L(P,P) \ RHS(P); or MAXITER = 20; TOL = 1e-13; RESTART = []; [ML,MU] = ilu(L(P,P),struct(’type’,’ilutp’,’droptol’,1e-6)); portion of system matrix after applying MATLAB symrcm u(P) = gmres(L(P,P),RHS(P),RESTART,TOL,MAXITER,ML,MU); 8 / 21
t+1D Advection Example solution in space-time domain PDE : u t = u x + F ( x , t ) ( x , t ) ∈ [ X ( t ) , 1) × (0 , T ] IC : u ( x , 0) = f ( x ) BC : u (1 , t ) = g ( t ) 1+( ε r ) 2 . r 2 = ( x − x i ) 2 + ( t − t i ) 2 1 √ IMQ-RBF: f ( x ) = e − 10( x − 0 . 15+0 . 35 y ) 2 Get RBF-QR diffmat from Elisabeth’s website. 10 -1 n loc = 50 , ε =3 0 10 10 -2 −1 10 10 -3 slope = -1.227 D t − D x −4.1 F −2 10 -4 10 ||.|| ∞ ||.|| = u 10 -5 −3 10 f 10 -6 0 I Test case 1 −4 10 Test case 2 g 10 -7 Test case 3 Test case 4 −5 10 -8 10 10 20 30 40 50 5 10 15 √ N 8 / 21
t+1D Advection with Variable Speed Example solution in space-time domain PDE : u t = a ( x , t ) u x + F ( x , t ) ( x , t ) ∈ [ X ( t ) , 1) × (0 , T ] IC : u ( x , 0) = f ( x ) BC : u (1 , t ) = g ( t ) a ( x , t ) = exp((1 + t )(1 + cos(3 x )) f ( x ) = e − 10( x − 0 . 15+0 . 35 y ) 2 Bribe Varun for PHS diffmat. 10 -1 n loc = 50 , ε =3 0 10 10 -2 −1 10 10 -3 slope = -1.456 D t − diag( a ) D x −4.1 F −2 10 -4 10 ||.|| ∞ ||.|| = u 10 -5 −3 10 f 10 -6 0 I Test case 1 −4 10 g Test case 2 10 -7 Test case 3 Test case 4 −5 10 -8 10 10 20 30 40 50 5 10 15 √ N 9 / 21
Analyzing Stability ? Let’s take a look at one step ( 2 levels) space-time global RBF method. ∆ x t + ∆ t . . . . . . x ′ x ′ x ′ x ′ x ′ 1 2 n − 1 n j t . . . . . . x 1 x 2 x j x n − 1 x n a b for a simple 1-D advection equation ∂ u ∂ t = ∂ u PDE: for x ∈ [ a , b ) ∂ x IC: u ( x , 0) = u 0 ( x ) when t = 0 BC: u ( b , t ) = g ( t ) at x = b n n � � λ ′ j φ ( ε � x − x ′ u ( x ) = λ j φ ( ε � x − x j � ) + j � ) , j =1 j =1 where { x j } and { x j } are centers at the old time level and new time level respectively. 10 / 21
∆ x t + ∆ t . . . . . . x ′ x ′ x ′ x ′ x ′ 1 2 n − 1 n j t . . . . . . x 1 x 2 x j x n − 1 x n a b Our goal is to find the unknowns { λ j } and { λ ′ j } . This can be done by enforcing initial and boundary data and satisfying the PDE at the interior points that lead to solving system of linear equations λ 1 u o ( x 1 ) A B . . . . . . λ n u o ( x n ) = λ ′ 0 1 C D . . . . . . λ ′ g ( t ) n 11 / 21
λ 1 u o ( x 1 ) A B . . . . . . λ n u o ( x n ) = λ ′ 0 1 C D . . . . . . λ ′ g ( t ) n The block matrices A , B , C , D are all n × n matrices with elements: A ij = φ ( ε � x i − x j � ) B ij = φ ( ε � x i − x ′ j � ) C ij = L φ ( ε � x ′ i − x j � ) D ij = L φ ( ε � x ′ i − x ′ j � ) for all i , j = 1 , · · · , n and L := ∂ ∂ t − ∂ ∂ x . The last row C and D must be slightly modified to satisfy the boundary condition at x ′ n = b . 12 / 21
Amplification Matrix and Stability Region The process of marching in time to the new time level is given by u ( x ′ 1 ) u ( x 1 ) . . . = . G . . u ( x ′ n ) u ( x n ) where A � � A � − 1 � G = � B � B I , C D 0 and I is an n × n identity matrix. The method is numerically stable if spectral radius ρ ( G ) < 1. −1 2 −2 1.5 IMQ, g ( t ) = 0, N = 50: to avoid blowing up the −3 −log( ε ) solution, the ratio of ∆ t / ∆ x vs shape parameter 1 ε must be away from the darker alley in the the −4 stability region, i.e we want to avoid ρ ( G ) ≥ 1 −5 0.5 −6 −4 −3 −2 −1 0 1 log( ∆ t / ∆ x) 13 / 21
Recommend
More recommend
