Solving the high-dimensional Vlasov equation with deal.II and - PowerPoint PPT Presentation

Solving the high-dimensional Vlasov equation with deal.II and hyper.deal Eighth deal.II Users and Developers Workshop Peter Munch 12 , Katharina Kormann 3 , Martin Kronbichler 1 1 Institute for Computational Mechanics, Technical University of Munich, Germany 2 Institute of Materials Research, Materials Mechanics, Helmholtz-Zentrum Geesthacht, Germany 3 Max Planck Institute for Plasma Physics and Department of Mathematics, Technical University of Munich, Germany May 26, 2020

Overview Vlasov equation: non-linear, high-dimensional, hyperbolic PDE � ∇ � � � � ∂ f ∂ f � v x ∂ t + � x f + � a ( t , f ,� x ,� v · ∇ � v ) · ∇ � v f = 0 ↔ ∂ t + · f = 0 � a ( t , f ,� x ,� v ) ∇ � v ... with additional term with derivative of the solution with respect to � v Table of contents: 1. Application: computational plasma physics 2. Discretization with discontinuous Galerkin methods 3. Solving with deal.II and hyper.deal 9.2 = new deal.II feature (release 9.2) 4. Extension to Vlasov–Poisson equation 1/26

Overview (cont.) hyper.deal ◮ FEM library with specialized algorithms for high dimensions (2 ≤ d ≤ 6) ◮ deal.II-based & 9.2 -ready ◮ freely available under LGPL 3.0 license ◮ hosted at https://github.com/hyperdeal/hyperdeal ◮ 2 tutorials: examples → advection & examples → vlasov poisson Munch, Kormann, Kronbichler, hyper.deal: An efficient, matrix-free finite-element library for high-dimensional partial differential equations, arXiv:2002.08110, 2020 2/26

Part 1: Application: computational plasma physics 3/26

Plasma physics Goal Describe the evolution of a plasma and its interaction in magnetic fields. A field of application is fusion energy research, in which the plasma in fusion reactors (e.g., tokamak and stellarator) are investigated. Mathematical descriptions: ⊲ n-body problem 1. particle model: description of the motion of each particle ∂ 2 x i ∂ t 2 = q i � � � v × � E ( t ,� x )+ � B ( t ,� x ) m i 2. kinetic model: described by a distribution function f , which evolves according to the Vlasov equation coupled to a system of Maxwell’s equations 3. fluid model: coupling of the Navier–Stokes equations to a system of Maxwell’s equations 4/26

Vlasov–Maxwell/Poisson equations Vlasov equation: with a single particle species with charge q and mass m ∂ f � � v ) = q � v × � ∂ t + � x f + � a ( t , f ,� x ,� � a ( t , f ,� x ,� E ( t ,� x )+ � B ( t ,� v · ∇ � v ) · ∇ � v f = 0 x ) m which is coupled to the Maxwell’s equations (or in simple cases to the Poisson equation) for the self-consistent fields. 5/26

Part 2: Discretization with discontinuous Galerkin methods 6/26

Discontinuous Galerkin discretization ◮ short-hand notation: � � � ∇ � � � ∂ f v x ∂ t + � with � a · ∇ f = 0 a := and ∇ := � a ∇ � v ◮ discontinuous Galerkin discretization expressed in reference space: � � � � � � g , |J | ∂ f J − T ∇ � af ) ∗ � g , |J | � � − ξ g , |J | + n · ( = 0 af ∂ t Ω ( e ) Ω ( e ) Γ ( e ) 0 0 0 7/26

Part 3: Solving with deal.II and hyper.deal 8/26

Limitation of with deal.II for high dimensions ⊲ step-67 To solve an advection equation in deal.II : 9.2 ◮ distributed triangulation: ◮ parallel::distributed::Triangulation<dim> ◮ parallel::fullydistributed::Triangulation<dim> 9.2 ◮ FE DGQ<dim>(k) and QGauss<dim>(k+1) ◮ MatrixFree<dim> infrastructure Problem: dim ≤ 3! In our case, dim=2, 3, 4 , 5 , ...! Question: What can we reuse from deal.II? Solution: tensor product of two (deal.II) triangulations! → hyper.deal 9/26

hyper.deal: triangulation T v T := T � x ⊗T � c v v c := ( c � x , c � v ) T x c x 10/26

hyper.deal: partitioning T v ◮ partition T x and T v independently c v ◮ checkerboard T x partitioning c x T f ( i , j ) := T i x ⊗T j f ( i , j ) := j · p � x + i with � � v 11/26

hyper.deal: partitioning (cont.) T v G 10 ◮ number of neighbors: L 10 larger c v ◮ neighboring cells: T x disjunct c x 12/26

hyper.deal: partitioning (cont.) MPI_COMM_WORLD column comm 24 25 column comm global comm 2 3 3 18 19 20 21 22 23 ◮ virtual topology 0 1 2 12 13 14 15 16 17 ◮ shared memory (MPI-3.0) row comm ◮ consensus algorithms 9.2 1 6 7 8 9 10 11 5 0 0 1 2 3 4 5 sm comm sm comm 2 3 4 0 1 2 3 4 5 0 1 row comm 13/26

hyper.deal: shape functions and quadrature ◮ due to tensor-product structure, an extension of the shape function to higher dimensions is simple: P d � x + d � = P d � k ⊗P d � = P 1 k ⊗···⊗P 1 ⊗P 1 k ⊗···⊗P 1 v x v k k k k � �� � � �� � × d � × d � x v ... with P d � and P d � k = P 1 k ⊗···⊗P 1 k = P 1 k ⊗···⊗P 1 x v k k � �� � � �� � × d � × d � x v ◮ the same is true for the quadrature rules ◮ as a consequence, the mapping data from lower-dimensional spaces (e.g., J x and J v ) from deal.II can be reused 14/26

hyper.deal: matrix-free infrastructure ◮ construct two dealii::MatrixFree instances ◮ loop over a pair of cell batches ◮ cell loop, face-centric loop, element-centric loop 9.2 ◮ access mapping information via FEEvaluation and FEFaceEvaluation 15/26

hyper.deal: matrix-free infrastructure (cont.) ∀ c := ( c � x , c � v ) hyperdeal::FEEvaluation reinit(cell) hyperdeal::MatrixFree read dof values(src) advection cell operation: ∀ q := ( q � x , q � v ) distribute local to global(dst) loop element centric() : void // ECL loop() : void // FCL get value < • > (val, q, qx, qv) f ( � u ) = J − 1 c , q |J q | w q � u , � u ∈ R d submit gradient < • > (val, q, qx, qv) 2 2 dealii::FEEvaluation dealii::MatrixFree � J − 1 � 0 c � x , c � u ∈ R d f ( � u ) = |J c � x | w q � x |J c v , q v | w q � v � u , � ShapeInfo x x , q � quadrature point(q) J − 1 0 v , q � MappingInfo c � v JxW(q) jacobian(q) 16/26

hyper.deal: matrix-free infrastructure (cont.) deal.II-like interface: // create hyperdeal::FEEvaluation FEEvaluation<dim_x, dim_v, degree, n_points, Number, VectorizedArrayType> phi( matrix_free, dof_no_x, dof_no_v, quad_no_x, quad_no_v); // configure underlying FEEvaluation // loop over cells (i.e. cell pairs) matrix_free.cell_loop([&](const auto &, auto & dst, const auto &src, const auto cell){ // reinit phi.reinit(cell); phi.read_dof_values(src); // evaluate // loop over quadrature points for (unsigned int qv = 0, q = 0; qv < phi.n_q_points_v; qv++) for (unsigned int qx = 0; qx < phi.n_q_points_x; qx++, q++) { // operation on quadrature point level const auto q_point = phi.get_quadrature_point(qx, qv); /*...*/ } // integrate phi.distribute_local_to_global(dst); }, dst, src); 17/26

hyper.deal: vectorization ◮ vectorization over a batch of elements ◮ constructing a batch of elements × VectoriedArray<T, 1> VectoriedArray<T, N> = VectoriedArray<T, N> � �� � � �� � x cell batch → vectorized v cell batch → not vectorized ... via the new template argument of VectorizedArray 9.2 ◮ setup of internal data structures of x-/v-MatrixFree appropriately MatrixFree<dim, T, VT> ... via the new template argument of VT := VectorizedArray<T, N> 9.2 18/26

hyper.deal: node-level performance Node-level performance (SuperMUC-NG) ◮ setup: Single operator application [GDoF/s] 4 tensor product of two hyper-cubes L2/L3 (see: examples → advection ) 2 ◮ challenge: temporal data size of cell batch O ( k d ) 1 k=3 + d=6 : 1.4GDoF/s ◮ outlook: vectorization within elements L1/L2 0.5 k = 2 k = 3 k = 4 k = 5 0.25 8 64 512 4096 32768 for different dimensions DoFs per cell 19/26

hyper.deal: strong and weak scaling Strong & weak scaling (SuperMUC-NG, d=6) 4.4TDoFs (Cartesian) 10 0 Time of 1 RK stage [s] 1.1GDoFs/node largest simulation: ◮ 4 . 4 · 10 12 = ( 2 . 1 · 10 6 ) 2 = 128 6 DoFs 268MDoFs/node 10 − 1 17.2GDoFs challenges: 2.1GDoFs ◮ communication pattern 10 − 2 ◮ communication data volume alternative partitioning 10 − 3 1 4 16 64 256 1024 3072 Nodes ( × 48 CPUs) 20/26

Part 4: Extension to Vlasov–Poisson equation 21/26

Vlasov-Poisson equation Vlasov equation for electrons in a neutralizing background in the absence of magnetic fields: � � ∂ f � v ∂ t + · ∇ f = 0 , − � E ( t ,� x ) where the electric field can be obtained from the Poisson problem: � − ∇ 2 � ρ ( t ,� x ) = 1 − f ( t ,� x ,� v ) d v , x φ ( t , x ) = ρ ( t ,� x ) , E ( t ,� x ) = − ∇ � x φ ( t ,� x ) . � ⊲ full code online: examples → vlasov poisson 22/26

Coupling with a deal.II-based Poisson solver ◮ in each Runge-Kutta step solve: � ρ ( t ,� f ( t ,� x ,� x ) = 1 − v ) d v � � ∇ 2 � x φ ( t , x ) = − ρ ( t ,� x ) ∂ f v � ∂ t + · ∇ f = 0 − � E ( t ,� x ) � E ( t ,� x φ ( t ,� x ) = − ∇ � x ) Ω x × Ω v → hyper.deal Ω x → deal.II → geometric multigrid ... the same approach is applicable to the Vlasov–Maxwell equations! 23/26

Part 5: Conclusions & outlook 24/26