Slappy: A BGK-Semi-Lagrangian parallel python solver on non-conforming patches

The Python library SLaPPy was developed to solve transport equa- tions on three-dimensional complex geometries. The geometries are divided into sub-domains called patches. SLaPPy is based on the backwards Semi-Lagrangian scheme. It uses local continuous Gauss- Lobatto interpolation inside the patches. The interpolation is dis- continuous between patches, thus it handles overlapping and non- conforming patches. The code is written in Python for the high-level data structures and OpenCL for the low-level computations and par- allelism.
The approximated BGK method allow us to split complex kinetic equations into independent transport equations, solved by the core of the SLaPPy code, and a relaxation step, which was added and paral- lelized in the code. Thus, we will present the results obtained with the library on linear and circular advections, isotropic diffusion and MHD models, on different geometries formed by non continuous patch.

Transcript

1. Application of the approximated BGK method on a Semi-Lagrangian parallel python solver on non-conforming patches Emmanuel Franck1, Philippe Helluy1 Laura S. Mendoza1 1Inria Grand-Est, TONUS TEAM NUMKIN 2018, Garching, Germany NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 1
2. Slappy: the objectives The motivation is to have a code capable of: solving transport equations § backwards Semi-Lagrangian method with local interpolation treat complex 3D geometries § domain decomposition on non conforming patches being a portable and fast code § Python for the high-level data structures and OpenCL for the low-level computations and parallelism. Application to complex models The Semi-Lagrangian Parallel Python solver on non-conforming patches NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 2
3. The Backward Semi-Lagrangian Method We consider the advection equation ∂f (x, t) ∂t + a(x, t) · ∇xf (x, t) = 0, f (x, 0) known (1) The scheme: Fixed grid on phase-space Method of characteristics : ODE −→ origin of characteristics Density f is conserved along the characteristics i.e. f n+1(xi) = f n(X(tn; xi, tn+1)) (2) Interpolate on the origin using known values of the previous time step at mesh points (initial distribution f 0 known). tn+1 t tn xi X(tn; xi , tn+1) NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 3
4. The Backward Semi-Lagrangian Method* We consider the advection equation ∂f (x, t) ∂t + a(x, t) · ∇xf (x, t) = 0, f (x, 0) known (3) The scheme: Fixed grid on space § Mesh refined in Nc cells and each cell L discretized into Nd points (Gauss-Lobatto quadrature) NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 4
5. The Backward Semi-Lagrangian Method* We consider the advection equation ∂f (x, t) ∂t + a(x, t) · ∇xf (x, t) = 0, f (x, 0) known (3) The scheme: Fixed grid on space § Mesh refined in Nc cells and each cell L discretized into Nd points (Gauss-Lobatto quadrature) Interpolation step § Local Lagrange interpolation in each cell L using the Lagrange polynomial basis function ϕ of degree d − 1 ˜ f (x, t) = Nd−1 j=0 fL,j (t)ϕL,j (x), x ∈ L NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 4
6. The Backward Semi-Lagrangian Method* We consider the advection equation ∂f (x, t) ∂t + a(x, t) · ∇xf (x, t) = 0, f (x, 0) known (3) The scheme: Fixed grid on space § Mesh refined in Nc cells and each cell L discretized into Nd points (Gauss-Lobatto quadrature) Interpolation step § Local Lagrange interpolation in each cell L using the Lagrange polynomial basis function ϕ of degree d − 1 ˜ f (x, t) = Nd−1 j=0 fL,j (t)ϕL,j (x), x ∈ L Advantages: Easy to have high-degree interpolations Possibility to have non conforming and/or overlapping patches Data needed for interpolation close-by memory-wise NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 4
7. Definition of patches The user defines the domain as a set of patches A patch is defined by a direct coordinate transformation τ, which can be: analytical § No need to store points coordinates obtained from pygmsh § Coordinates needed in a rough refinement Discretization done with a structured grid Interpolation step done in the referential domain § The inverse function is either given analytically or computed with a Newton τ NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 5
8. Structure of the code Slappy is written in Python for the high-level data structures and OpenCL for the low-level computations and parallelism. Python Object-oriented code structure Mesh creation (analytical or Pygmsh) Time and patches loops Creation and distribution of OpenCL kernels and buffers Post evaluation (vtk, ...) OpenCL (C) Kernels written in C99 Execution on heterogeneous platforms (CPU/GPU) Paralellized computation on the Gaussian points Handles device to device memory transfers Interfaced with PyOpenCL NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 6
9. Advection equation on a 3D non-conforming mesh ∂f (x, t) ∂t + a · ∇xf (x, t) = 0, f (x, 0) known, with a = (−0.2, −0.2, 0). Refinement 20, 15, 20, and degree 3 for each patch. NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 7
10. Advection equation on a 3D non-conforming mesh ∂f (x, t) ∂t + a · ∇xf (x, t) = 0, f (x, 0) known, with a = (−0.2, −0.2, 0). Simulation parameters t = 0.0, ∆t = 0.2, tmax = 5 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 8
11. Advection equation on a 3D non-conforming mesh ∂f (x, t) ∂t + a · ∇xf (x, t) = 0, f (x, 0) known, with a = (−0.2, −0.2, 0). Simulation parameters t = 1., ∆t = 0.2, tmax = 5 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 9
12. Advection equation on a 3D non-conforming mesh ∂f (x, t) ∂t + a · ∇xf (x, t) = 0, f (x, 0) known, with a = (−0.2, −0.2, 0). Simulation parameters t = 1.8, ∆t = 0.2, tmax = 5 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 10
13. Advection equation on a 3D non-conforming mesh ∂f (x, t) ∂t + a · ∇xf (x, t) = 0, f (x, 0) known, with a = (−0.2, −0.2, 0). Simulation parameters t = 2.6, ∆t = 0.2, tmax = 5 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 11
14. Advection equation on a 3D non-conforming mesh ∂f (x, t) ∂t + a · ∇xf (x, t) = 0, f (x, 0) known, with a = (−0.2, −0.2, 0). Simulation parameters t = 5, ∆t = 0.2, tmax = 5 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 12
15. Linear Advection - Study of convergence Deg Raf min(f) max(fn)/max(f0) L∞ L2 1 10 4e-07 14.17 2.32e-06 1.43e-05 1 20 4e-12 3.7 3.9e-08 4.49-07 1 40 8e-17 1.68 3.51e-10 1.01e-08 2 10 -e-04 1.6 2.46e-06 1.38e-05 2 20 -0.01 1.08 1.04e-08 1.08e-07 2 40 -2e-04 1.01 1.65e-11 4.82e-10 3 10 -0.37 1.11 1.81e-07 1.93e-06 3 20 -6e-05 1.005 1.57e-10 4.10e-09 3 40 4.33e-12 1.006 1.39e-13 8.61e-12 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 13
16. Linear Advection - Scalability Scalability study on CPU: #Units #Points Time Ite. Exec. Time (s) 1 921 600 25 152 2 921 600 25 71 4 921 600 25 37 8 921 600 25 25 8 1 843 200 25 45 16 1 843 200 25 25 GPU (NVIDIA - CUDA, Tesla K80) vs CPU (64 cores): Device #Points Time Ite. Exec. Time (s) GPU 17 203 200 25 75 CPU 17 203 200 25 225 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 14
17. Generalization - Approximated BGK model Let a hyperbolic model ∂tU + ∂xF(U) = 0. (4) Then, the approximated BGK model ∂tf + Λ∂xf = 1 ε (f eq(U) − f ) with Λ matrix of velocities and with the condition Pf eq(U) = U, PΛf eq(U) = F(U). tends to (4) when ε → 01. 1S. Jin and Z. Xin. Communications on pure and applied mathematics 48.3 (1995), pp. 235–276. NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 15
18. Approximated BGK model - Discretization Let us discretize the approximated BGK model ∂tf + Λ∂xf = 1 ε (f eq(U) − f ). (5) Then, the splitting scheme T∆t ◦ R∆t with Advection T∆t : f n+1 = f n(x − Λ∆t) Relaxation R∆t : f n+1 = f n + ω(f eq(U) − f n) with ω = ∆t θ∆t + ε is consistent with (5). NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 16
19. Approximated BGK model - Discretization Let us discretize the approximated BGK model ∂tf + Λ∂xf = 1 ε (f eq(U) − f ). (5) Then, the splitting scheme T∆t ◦ R∆t with Advection T∆t : f n+1 = f n(x − Λ∆t) Relaxation R∆t : f n+1 = f n + ω(f eq(U) − f n) with ω = ∆t θ∆t + ε is consistent with (5). In Slappy: existent advection module, added (parallelized) relaxation NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 16
20. Numerical results - Diffusion Model of isotropic diffusion ∂tU − ∂xxU − ∂yyU = 0 (6) D2Q4 - Approximated BGK model § Velocity set : Λ = {(1, 0)T , (0, 1)T , (−1, 0)T , (0, −1)T } § ∀i ∈ 1, · · · , 4, f eq i (U) = 1 4 U § U = 4 i=0 fi NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 17
21. Numerical results - Diffusion Simulation parameters : § 5 patch disk, refinement (15, 15) for external patches and (30, 30) for internal patch § Interpolation basis function degree 3 § tmax = 0.25, ∆t = 0.01 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 18
22. Numerical results - Diffusion Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 0.02, ∆t = 0.001 § t = 0 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 19
23. Numerical results - Diffusion Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 0.02, ∆t = 0.001 § t = 0.002 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 20
24. Numerical results - Diffusion Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 0.02, ∆t = 0.001 § t = 0.005 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 21
25. Numerical results - Diffusion Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 0.02, ∆t = 0.001 § t = 0.01 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 22
26. Numerical results - Diffusion Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 0.02, ∆t = 0.001 § t = 0.02 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 23
27. Numerical results - Circular Advection Model of circular advection ∂tU + ∇xA(x)U = 0 (7) D2Q4 - Approximated BGK model § Velocity set : Λ = {(1, 0)T , (0, 1)T , (−1, 0)T , (0, −1)T } § ∀i ∈ 1, · · · , 4, f eq i (U) = 1 4 U + 1 2λ2 i (Af × λi) § U = 4 i=0 fi NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 24
28. Numerical results - Circular Advection Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 1.5, ∆t = 0.001 § t = 0.25 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 25
29. Numerical results - Circular Advection Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 1.5, ∆t = 0.001 § t = 0.3 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 26
30. Numerical results - Circular Advection Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 1.5, ∆t = 0.001 § t = 0.6 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 27
31. Numerical results - Circular Advection Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 1.5, ∆t = 0.001 § t = 0.9 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 28
32. Numerical results - Circular Advection Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 1.5, ∆t = 0.001 § t = 1.2 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 29
33. Numerical results - Circular Advection Simulation parameters : § refinement [(15, 15)4, (30, 30)] § Interpolation basis function degree 3 § tmax = 1.5, ∆t = 0.001 § t = 1.5 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 30
34. Numerical results - Isothermal Euler Model of isothermal Euler ∂tρ + ∇ · (ρu) = 0 ∂tρu + ∇ · (ρu × u + c2ρId) = 0 § Refinement [(15, 15)4, (20, 20)] § Interpolation basis function degree 3 § tmax = 0.5, ∆t = 0.0005 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 31
35. Numerical results - Isothermal Euler § t = 0 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 32
36. Numerical results - Isothermal Euler § t = 0.O625 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 33
37. Numerical results - Isothermal Euler § t = 0.125 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 34
38. Numerical results - Isothermal Euler § t = 0.1875 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 35
39. Numerical results - Isothermal Euler § t = 0.25 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 36
40. Numerical results - Isothermal Euler § t = 0.3125 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 37
41. Numerical results - Isothermal Euler § t = 0.375 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 38
42. Numerical results - Isothermal Euler § t = 0.375 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 39
43. Numerical results - Isothermal Euler § t = 0.4375 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 40
44. Numerical results - Isothermal Euler § t = 0.5 NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 41
45. Slappy: the objectives Slappy is capable of: solving transport equations § backwards Semi-Lagrangian method with local interpolation treat complex 3D geometries § domain decomposition on non conforming patches being a portable and fast code § Python for the high-level data structures and OpenCL for the low-level computations and parallelism. Application to complex models Perspectives: § Toroidal geometry (on going) § MHD, Non isotropic diffusion, ... § Hybrid parallelization NUMKIN 2018 – Laura S. Mendoza (laura.mendoza@inria.fr) Wednesday 24th October, 2018 42