Preconditioning and nonlinear time solvers for the JOREK MHD code - PowerPoint PPT Presentation

Physical context and models JOREK code and time solvers Preconditioning Preconditioning and nonlinear time solvers for the JOREK MHD code E. Franck, A. Lessig, M. H¨ olzl, E. Sonnendr¨ ucker Max Planck Institute of Plasma physics, Garching, Germany 10th AIMS Conference, special session 121, Madrid 10 July 2014 E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Outline 1 Physical context and models 2 JOREK code and time solvers 3 Preconditioning E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Physical context and models E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Magnetic Confinement Fusion Fusion DT : At sufficiently high energies, deuterium and tritium can fuse to Helium. A neutron and 17.6 MeV of free energy are released. At those energies, the atoms are ionized forming a plasma. E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Magnetic Confinement Fusion Fusion DT : At sufficiently high energies, deuterium and tritium can fuse to Helium. A neutron and 17.6 MeV of free energy are released. At those energies, the atoms are ionized forming a plasma. Magnetic confinement : The charged plasma particles can be confined in a toroidal magnetic field configuration, for instance a tokamak. Figure : Tokamak E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Plasma instabilities Edge localized modes (ELMs) are periodic instabilities occurring at the edge of tokamak plasmas. They are associated with strong heat and particle losses which could damage wall components in ITER by large heat loads. Aim : Detailed non-linear modeling and simuation (MHD models) can help to understand and control ELMs better. Final Density Initial Density E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning MHD model The full resistive MHD model is given by ∂ t ρ + ∇ · ( ρ v ) = ∇ · ( D ∇ ρ ) + S p       ρ∂ t v + ρ v · ∇ v + ∇ P = J × B + ν △ v        ∂ t P + v · ∇ P + γ P ∇ v = ∇ · ( K ∇ T ) + S h     ∂ t B = −∇ × E = ∇ × ( v × B ) − η ∇ × J         ∇ · B = 0 Magnetic quantities : B the magnetic field, E the electric field and J = ∇ × B the current. Hydrodynamic quantities : ρ the density, v the velocity, T the temperature, and P = ρ T the pressure. The terms K and D are anisotropic diffusion tensors. Source terms: S h is a heat source, S p is a particle source. E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Reduced MHD: assumptions and principle of derivation Aim : Reduce the number of variables and eliminate the fast magnetosonic waves. We consider the cylindrical coordinate ( R , Z , φ ) ∈ Ω × [0 , 2 π ] Reduced MHD: Assumptions B = F 0 R e φ + 1 R ∇ ψ × e φ v = − R ∇ u × e φ + v || B with u the electrical potential, ψ the magnetic poloidal flux, v || the parallel velocity. To avoid high order operators we introduce the vorticity w = △ pol u and the toroidal current j = △ ∗ ψ = R 2 ∇ · ( 1 R 2 ∇ pol ψ ). Derivation: we plug B and v in the equations + some computations. For the equations on u and v || we use the following projections e φ · ∇ × R 2 ( ρ∂ t v + ρ v · ∇ v + ∇ P = J × B + ν △ v ) and B · ( ρ∂ t v + ρ v · ∇ v + ∇ P = J × B + ν △ v ) . E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Reduced MHD without v || : simple model Example of model: case where v || = 0. ∂ t ψ = R [ ψ, u ] − F 0 ∂ φ u + η ( T )( j + 1  R 2 ∂ φφ ψ )         ρ ∇ pol ( ∂ t u )) = 1 ρ w , u ] + [ ψ, j ] − F 0  ρ ] + [ R 2 ˆ  2 [ R 2 ||∇ pol u || 2 , ˆ R ∂ φ j − [ R 2 , P ] R ∇ · (ˆ       + ν R ∇ · ( ∇ pol w )        R 2 j − ∇ · ( 1 1 R 2 ∇ pol ψ ) = 0        w − ∇ · ( ∇ pol u ) = 0         ∂ t ρ = R [ ρ, u ] + 2 ρ∂ Z u + ∇ · ( D ∇ ρ )         ∂ t T = R [ T , u ] + 2( γ − 1) T ∂ Z u + ∇ · ( K ∇ T ) ρ = R 2 ρ . with ˆ D and K are anisotropic diffusion tensors (in the direction parallel to B ). η ( T ) is the physical resistivity. ν is the viscosity. E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Main result: energy estimate Correct reduced model : estimation on the energy conservation or dissipation. Model with parallel velocity: We assume that the boundary conditions are correctly chosen. The fields are defined by B = F 0 R e φ + 1 R ∇ ψ × e φ and v = − R ∇ u × e φ + v || B . For the model associated with these fields we obtain η |△ ∗ ψ | 2 η |∇ pol ( ∂ φ ψ d � � � � R 2 ) | 2 − ν |△ pol u | 2 E ( t ) = − − R 2 dt Ω Ω Ω Ω with E ( t ) = | B | 2 + ρ | v | 2 1 + γ − 1 P the total energy. 2 2 The implemented models approximately conserve energy. For exact energy conservation, some neglected cross-terms between poloidal and parallel velocity have to be added which might be important in the non-linear phase. Theoretical and numerical stability for the reduced MHD models in JOREK code , E. Franck, M. H¨ olzl, A. Lessig, E. Sonnendr¨ ucker, in redaction E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Jorek code and time solvers E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Description of the JOREK code I JOREK: Fortran 90 code, parallel (MPI+OpenMP) + algebraic libraries (Pastix, MUMPS ...) Initialization Determine the equilibrium Define the boundary of the computational domain Create a first grid which is used to compute the aligned grid Compute ψ ( R , Z ) in the new grid. Compute equilibrium Solve the Grad-Shafranov equation � 1 + ∂ 2 ψ R ∂ ∂ψ � ∂ Z 2 = − R 2 ∂ p ∂ψ − F ∂ F Figure : unaligned grid ∂ R R ∂ R ∂ψ E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Description of the JOREK code II Computation of aligned grid Identification of the magnetic flux surfaces Create the aligned grid (with X-point) Interpolate ψ ( R , Z ) in the new grid. Recompute equilibrium of the new grid. Perturbation of the equilibrium (small perturbations of non principal harmonics). Time-stepping (full implicit) Poloidal discretization : 2D Cubic Bezier finite elements. Toroidal discretization : Fourier expansion. Construction of the matrix and some profiles (diffusion tensors, sources terms). Solve linear system. Figure : Aligned grid Update solutions. E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Time scheme in JOREK code The model is ∂ t A ( U ) = B ( U , t ) For time stepping we use a Crank Nicholson or Gear scheme : (1 + ζ ) A ( U n +1 ) − ζ A ( U n ) + ζ A ( U n − 1 ) = θ ∆ tB ( U n +1 ) + (1 − θ )∆ tB ( U n ) Defining G ( U ) = (1 + ζ ) A ( U ) − θ ∆ tB ( U ) and b ( U n , U n − 1 ) = (1 + 2 ζ ) A ( U n ) − ζ A ( U n − 1 ) + (1 − θ )∆ tB ( U n ) we obtain the nonlinear problem G ( U n +1 ) = b ( U n , U n − 1 ) First order linearization � ∂ G ( U n ) � δ U n = − G ( U n ) + b ( U n , U n − 1 ) = R ( U n ) ∂ U n with δ U n = U n +1 − U n , and J n = ∂ G ( U n ) the Jacobian matrix of G ( U n ). ∂ U n E. Franck and al. Nonlinear time solvers for Jorek MHD code

Physical context and models JOREK code and time solvers Preconditioning Linear Solvers Linear solver in JOREK: Left Preconditioning + GMRES iterative solver. Principle of the preconditioning step: Replace the problem J k δ U k = R ( U n ) by P k ( P − 1 J k ) δ U k = R ( U n ). k Solve the new system with two steps P k δ U ∗ k = R ( U n ) and ( P − 1 J k ) δ U k = δ U ∗ k k If P k is easier to invert than J k and P k ≈ J k the linear solving step is more robust and efficient. Construction and inversion of P k P k : diagonal block matrix where the sub-matrices are associated with each toroidal harmonic. Inversion of P k : We use a LU factorization and invert exactly each subsystem. This preconditioning is based on the assumption that the coupling between the toroidal harmonics is weak. In practice for some test cases this coupling is strong in the nonlinear phase. E. Franck and al. Nonlinear time solvers for Jorek MHD code