Atmospheric Dynamics with Polyharmonic Spline RBFs Greg Barnett - PowerPoint PPT Presentation

Photos placed in horizontal position with even amount of white space between photos and header Atmospheric Dynamics with Polyharmonic Spline RBFs Greg Barnett Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administrati on under contract DE-NA0003525. SAND2017-8166 C

Outline ▪ Polyharmonic Spline (PHS) RBFs with Polynomials ▪ 1D Example ▪ Interpolation and Differentiation Weights ▪ Interpolation/Weights in 2D ▪ Weights on the Sphere ▪ Semi-Lagrangian Transport ▪ Governing Equations ▪ Limiter/Fixer ▪ Test Cases ▪ Results ▪ Eulerian Shallow Water Model ▪ Governing Equations ▪ Test Cases ▪ Results ▪ Conclusions and Future Work 2

Table of RBFs 𝜚 𝒚 , 𝒚 ∈ ℝ 𝑒 , 𝜁 ∈ ℝ Name (acronym) Multiquadric (MQ) 2 1 + 𝜁 𝒚 1 Inverse Quadratic (IQ) 2 1 + 𝜁 𝒚 1 Inverse Multiquadric (IMQ) 2 1 + 𝜁 𝒚 2 𝑓 − 𝜁 𝒚 Gaussian (GA) 𝒚 2𝑙+1 Polyharmonic Spline (PHS) 𝒚 2𝑙 log 𝒚 , 𝑙 ∈ ℕ ⋅ = ⋅ 2 3

Some RBFs in 1D Polyharmonic Spline RBFs Infinitely Differentiable RBFs 4

Example: Equi-spaced Interpolation in 1D PHS Basis Functions Approximation 5

Structure of the Linear System 2 1 𝑦 1 𝑦 1 𝑔 1 2 1 𝑦 2 𝑦 2 𝜈 1 𝑔 2 Least Squares Parabola: 2 𝜈 2 1 𝑦 3 𝑦 3 = 𝑔 3 𝜈 1 + 𝜈 2 𝑦 + 𝜈 3 𝑦 2 𝜈 3 2 𝑔 1 𝑦 4 𝑦 4 4 𝑔 2 1 𝑦 5 𝑦 5 5 2 3 4 1 𝑦 1 𝑦 1 𝑦 1 𝑦 1 𝜈 1 𝑔 1 2 3 4 1 𝑦 2 𝑦 2 𝑦 2 𝑦 2 𝜈 2 𝑔 2 Polynomial Interpolant: 2 3 4 𝜈 3 1 𝑦 3 𝑦 3 𝑦 3 𝑦 3 = 𝑔 𝜈 1 + 𝜈 2 𝑦 + 𝜈 3 𝑦 2 + 𝜈 4 𝑦 3 + 𝜈 5 𝑦 4 3 𝜈 4 2 3 4 𝑔 1 𝑦 4 𝑦 4 𝑦 4 𝑦 4 4 𝜈 5 𝑔 2 3 4 1 𝑦 5 𝑦 5 𝑦 5 𝑦 5 5 𝑦 1 − 𝑦 2 3 𝑦 1 − 𝑦 3 3 𝑦 1 − 𝑦 4 3 𝑦 1 − 𝑦 5 3 2 0 1 𝑦 1 𝑦 1 𝜇 1 𝑔 1 𝑦 2 − 𝑦 1 3 𝑦 2 − 𝑦 3 3 𝑦 2 − 𝑦 4 3 𝑦 2 − 𝑦 5 3 2 0 1 𝑦 2 𝑦 2 𝜇 2 𝑔 2 𝑦 3 − 𝑦 1 3 𝑦 3 − 𝑦 2 3 𝑦 3 − 𝑦 4 3 𝑦 3 − 𝑦 5 3 2 0 1 𝑦 3 𝑦 3 𝜇 3 𝑔 PHS Interpolant: 3 𝑦 4 − 𝑦 1 3 𝑦 4 − 𝑦 2 3 𝑦 4 − 𝑦 3 3 𝑦 4 − 𝑦 5 3 2 𝜇 4 0 1 𝑦 4 𝑦 4 𝑔 5 4 3 + 𝜈 1 + 𝜈 2 𝑦 + 𝜈 3 𝑦 2 𝜇 5 = 2 𝑦 5 − 𝑦 1 3 𝑦 5 − 𝑦 2 3 𝑦 5 − 𝑦 3 3 𝑦 5 − 𝑦 4 3 𝑔 0 1 𝑦 5 𝑦 5 5 ෍ 𝜇 𝑘 𝑦 − 𝑦 𝑘 𝜈 1 𝑘=1 0 1 1 1 1 1 0 0 0 𝜈 2 0 𝑦 1 𝑦 2 𝑦 3 𝑦 4 𝑦 5 0 0 0 𝜈 3 0 2 2 2 2 2 𝑦 1 𝑦 2 𝑦 3 𝑦 4 𝑦 5 0 0 0 6

Properties of Polyharmonic Splines ▪ PHS basis includes both RBFs and polynomials ▪ RBFs improve performance and allow the use of irregular nodes ▪ polynomials give convergence to smooth solutions (no saturation error) ▪ The interpolation problem is guaranteed to have a unique solution provided that polynomials are included up to the required degree, and the nodes are unisolvent. For 𝑙 = 1,2,3, … 𝒚 2𝑙 log 𝒚 ▪ 𝜚 𝒚 = ▪ Polynomials up to degree 𝑙 or higher ▪ 𝒚 2𝑙+1 𝜚 𝒚 = ▪ Polynomials up to degree 𝑙 or higher ▪ Rule of thumb for modest polynomial degrees: Twice as many RBFs as polynomials ▪ Condition number of PHS 𝐵 -matrix is invariant under rotation, translation, and uniform scaling ▪ No need to search for optimal shape parameter 7

Interpolation in 2D 𝒚 = 𝑦, 𝑧 𝑜 𝑜 𝑦 𝑗 , 𝑧 𝑗 and corresponding function values 𝑔 Given nodes , find a linear combination of RBF and 𝑗 𝑗=1 𝑗=1 polynomial basis functions that matches the data exactly. 1. Assume the appropriate form of the underlying approximation: 𝑜 𝑛 Φ 𝑦, 𝑧 = ෍ 𝜇 𝑘 𝜚 𝑘 𝑦, 𝑧 + ෍ 𝜈 𝑙 𝑞 𝑙 𝑦, 𝑧 , 𝑘=1 𝑙=1 where 𝜚 𝑘 𝑦, 𝑧 = 𝜚 𝑦 − 𝑦 𝑘 , 𝑧 − 𝑧 𝑘 . 2. Require Φ to match the data at each node: 𝑜 𝑛 Φ 𝑦 𝑗 , 𝑧 𝑗 = ෍ 𝜇 𝑘 𝜚 𝑘 𝑦 𝑗 , 𝑧 𝑗 + ෍ 𝜈 𝑙 𝑞 𝑙 𝑦 𝑗 , 𝑧 𝑗 = 𝑔 𝑗 , 𝑗 = 1,2,3, … , 𝑜. 𝑘=1 𝑙=1 3. Enforce regularity conditions on the coefficients 𝜇 𝑘 : 𝑜 𝜇 𝑘 𝑞 𝑙 𝑦 𝑘 , 𝑧 𝑘 = 0, 𝑙 = 1,2,3, … , 𝑛. ෍ 𝑘=1 4. Solve the symmetric linear system for 𝜇 𝑘 and 𝜈 𝑙 . 8

Differentiation Weights in 2D Interpolation Problem: 𝝁 𝐁 𝐐 𝝂 = 𝒈 𝑷 , 𝐐 𝑈 𝐏 𝑏 𝑗𝑘 = 𝜚 𝑘 𝑦 𝑗 , 𝑧 𝑗 = 𝜚 𝑦 𝑗 − 𝑦 𝑘 , 𝑧 𝑗 − 𝑧 𝑘 , 𝑗, 𝑘 = 1,2,3, … , 𝑜, 𝑞 𝑗𝑙 = 𝑞 𝑙 𝑦 𝑗 , 𝑧 𝑗 , 𝑗 = 1,2,3, … , 𝑜, 𝑙 = 1,2,3, … , 𝑛. Use 𝑀Φ 𝑦, ෤ ෤ 𝑧 to approximate 𝑀𝑔 𝑦, ෤ ෤ 𝑧 : 𝑜 𝑛 𝑀Φ 𝑦, ෤ ෤ 𝑧 = ෍ 𝜇 𝑘 𝑀𝜚 𝑘 𝑦, ෤ ෤ 𝑧 + ෍ 𝜈 𝑙 𝑀𝑞 𝑙 𝑦, ෤ ෤ 𝑧 𝑘=1 𝑙=1 −1 𝝁 𝐁 𝐐 𝒈 = 𝒄 𝝂 = 𝑷 , 𝒅 𝒄 𝒅 𝐐 𝑈 𝐏 weights 𝑐 𝑘 = 𝑀𝜚 𝑘 𝑦, ෤ ෤ 𝑧 , 𝑘 = 1,2,3, … , 𝑜, 𝑑 𝑙 = 𝑀𝑞 𝑙 𝑦, ෤ ෤ 𝑧 , 𝑙 = 1,2,3, … , 𝑛. 9

Derivative Approximation on the Sphere ▪ Given: nodes 𝒚, 𝒛, 𝒜 on the sphere and function values 𝒈 𝜖 𝜖 𝜖 ▪ Find: differentiation matrices (DMs) 𝐗 𝑦 , 𝐗 𝑧 , 𝐗 𝑨 to approximate 𝜖𝑦 , 𝜖𝑧 , 𝜖𝑨 ▪ Method: Use the fact that 𝛼𝑔 ෤ 𝑦, ෤ 𝑧, ǁ 𝑨 is tangent to the sphere at ෤ 𝑦, ෤ 𝑧, ǁ 𝑨 ▪ For each node, get orthogonal unit vectors 𝒇 ො 𝑦 and 𝒇 ො 𝑧 tangent to the sphere 𝜖 𝜖 ▪ Use 2D method to get matrices 𝐗 ො 𝑦 and 𝐗 ො 𝑧 that approximate 𝑦 and 𝜖 ො 𝜖 ො 𝑧 𝜖𝑔 𝜖𝑔 ▪ 𝛼𝑔 = 𝒇 ො 𝑦 + 𝒇 ො 𝑦 𝑧 𝜖 ො 𝜖 ො 𝑧 𝜖𝑔 𝜖𝑔 𝜖𝑔 𝜖 𝜖 ▪ 𝜖𝑦 = 𝛼𝑔 1 = 𝒇 ො 𝑦 + 𝒇 ො 𝑧 = 𝒇 ො 𝑦 + 𝒇 ො 𝑧 𝑔 𝑦 1 𝑧 1 𝑦 1 𝑧 1 𝜖 ො 𝜖 ො 𝜖 ො 𝜖 ො ▪ 𝐗 𝑦 = diag 𝒇 ො 𝑦 1 𝐗 ො 𝑦 + diag 𝒇 ො 𝑧 1 𝐗 ො 𝑧 ▪ 𝐗 𝑧 = diag 𝒇 ො 𝑦 2 𝐗 ො 𝑦 + diag 𝒇 ො 𝑧 2 𝐗 ො 𝑧 ▪ 𝐗 𝑨 = diag 𝒇 ො 𝑦 3 𝐗 ො 𝑦 + diag 𝒇 ො 𝑧 3 𝐗 ො 𝑧 10

Semi-Lagrangian Transport ▪ Governing Equations (velocity 𝒗 is a known function) 𝜖𝜍 ▪ 𝜖𝑢 = −𝒗 ⋅ 𝛼𝜍 − 𝜍𝛼 ⋅ 𝒗, (Eulerian, short time-steps) 𝐸𝑟 𝜖 ▪ 𝐸𝑢 = 𝜖𝑢 + 𝒗 ⋅ 𝛼 𝑟 = 0. (Semi-Lagrangian, long time-steps) ▪ Quasi-Monotone Limiter for 𝑟 (const. along flow trajectories) (𝑜) and 𝑁 𝑙 = max (𝑜) ▪ 𝑛 𝑙 = min 𝑟 𝑙 ℓ 𝑟 𝑙 ℓ ℓ ℓ (𝑜+1) = min 𝑟 𝑙 (𝑜+1) , 𝑁 𝑙 ▪ Set 𝑟 𝑙 (𝑜+1) = max 𝑟 𝑙 (𝑜+1) , 𝑛 𝑙 ▪ Set 𝑟 𝑙 𝑂 ▪ Mass Fixer tracerMass ≡ σ 𝑙=1 𝜍 𝑙 𝑟 𝑙 𝑊 𝑙 ▪ If tracerMass < initialMass, add mass in cells with 𝑟 𝑙 < 𝑁 𝑙 ▪ If tracerMass > initialMass, subtract mass from cells with 𝑟 𝑙 > 𝑛 𝑙 11

Semi-Lagrangian Transport 𝑢 𝑜+1 𝑢 𝑜 12

Pre-Processing ▪ Get DMs 𝐗 𝑦 , 𝐗 𝑧 , and 𝐗 𝑨 for the continuity equation ▪ Find index of the 𝑜 nearest neighbors to each node ▪ For each node 𝒚 𝑙 = 𝑦 𝑙 , 𝑧 𝑙 , 𝑨 𝑙 , write its 𝑜 neighbors in terms of two orthogonal unit vectors 𝒇 ො 𝑦 , 𝒇 ො 𝑧 tangent to the sphere at 𝒚 𝑙 , and one unit vector 𝒇 Ƹ 𝑨 normal to the sphere at 𝒚 𝑙 , so that 𝒚 𝑙 ℓ = ො 𝑦 𝑙 ℓ 𝒇 ො 𝑦 + ො 𝑧 𝑙 ℓ 𝒇 ො 𝑧 + Ƹ 𝑨 𝑙 ℓ 𝒇 Ƹ 𝑨 , ℓ = 1,2,3, … , 𝑜. −1 ▪ Set 𝐃 𝑙 = 𝐁 𝑙 𝐐 𝑙 𝐉 𝑜 𝐏 6×𝑜 , where 𝑈 𝐐 𝑙 𝐏 6×6 𝐁 𝑙 𝑗𝑘 = 𝜚 𝑦 𝑙 𝑗 − ො ො 𝑦 𝑙 𝑘 , ො 𝑧 𝑙 𝑗 − ො 𝑧 𝑙 𝑘 , 𝑗, 𝑘 = 1,2,3, … , 𝑜, 2 . 2 , ෝ 𝐐 𝑙 = 𝟐, ෝ 𝒚 𝑙 , ෝ 𝒛 𝑙 , ෝ 𝒚 𝑙 ෝ 𝒛 𝑙 , ෝ 𝒚 𝑙 𝒛 𝑙 13

Time Stepping 𝜖𝜍 𝜖𝑢 = −𝒗 ⋅ 𝛼𝜍 − 𝜍𝛼 ⋅ 𝒗 from 𝑢 𝑜 to 𝑢 𝑜+1 using several 1. Step explicit Eulerian time steps (RK3) 2. Step 𝒚 ′ = −𝒗 from 𝑢 𝑜+1 to 𝑢 𝑜 to get departure points (RK4) 3. Find the nearest fixed neighbor to each departure point 4. Use the corresponding pre-calculated cardinal coefficients 𝐃 𝑙 and the newly formed row-vectors 𝒄 𝑙 to get rows of the interpolation matrix 𝐗 𝐗 𝑙⋅ = 𝒄 𝑙 𝐃 𝑙 5. Update 𝑟 on fixed nodes using weights 𝐗 6. Cycle quasi-monotone limiter and mass-fixer until the tracer mass is nearly equal to the initial tracer mass (diff<1e-13) 7. Repeat 14