# Presentation: 938869

## Scalable Solvers in Petascale Electromagnetic Simulation

***Lie-Quan (Rich) Lee, ***Volkan Akcelik,
- Ernesto Prudencio, Lixin Ge
- Stanford Linear Accelerator Center
- Xiaoye Li, Esmond Ng
- Lawrence Berkeley National Laboratory

- Work supported by  DOE ASCR, BES & HEP Divisions under contract DE-AC02-76SF00515

## Overview

- Shape Determination/Optimization
- V. Akcelik, L. Lee (SLAC)‏
- T. Tautges, P. Knupp,  L. Diachin (ITAPS)‏
- O. Ghattas, E. Ng, D. Keyes (TOPS)‏
- Linear and Nonlinear Eigensolvers
- L. Lee(SLAC), X. Li, E. Ng, C. Yang (LBNL/TOPS)‏
- Scalable Linear Solvers
- L. Lee (SLAC), X. Li, E. Ng (TOPS)‏

## Shape Determination and Optimization

## Shape Determination and Optimization For SCRF Cavities

- Shape changes due to
- Fabrication errors
- Addition of stiffening rings
- Tuning for accelerating mode
- Change HOM Damping -> Beam quality

- Ring in the middle

- HOM Damping changes

- Tuning

## Least-squares Minimization

- Unknowns are shape deviation parameters
- Gauss-Newton with truncated-SVD
- Indefinite linear systems from KKT (deferred)‏

- Its forward problem
- is Maxwell eigenvalue
- problem

## Example 1 for ILC TDR Cavity

- Create a synthetic example, artificially deform a 3D 9 cell ILC cavity.
- Choose a set of parameters defining shape variations, in total 26 independent inversion parameters.
- Cell radius dr (x9) an cell length dz (x9)‏
- Iris radius (x8)‏
- Assign random values to these variables, and deform the cavity.
- Solve the Maxwell eigenvalue problem.
- Use the first 45 nonzero frequencies, and first 9 modes field distribution as the targeted values

## Results for Example 1

- The nonlinear solver converges within a handful of iterations
- Frequencies and Fields match remarkably
- Objective function decreases by 10e6
- The “target” and “inverted” cavity shapes are very close to each other

## Determining TDR Shape with Measured Frequencies 

- Experimental data for manufactured baseline ILC cavities from DESY
- The first 45 mode frequencies, and the first 9 monopole mode field distribution along the cavity axis
- 82 parameters: cell radius, length, tuning, warping, and iris radius

- Cell length error

- Cell radius error

- Deformed surface

- Elliptical shape

## Results

- Difference of Frequencies and Field values
- Red: inverted cavity - measured values
- Black/blue: ideal shape - measured values

- An article has been accepted by JCP

- MHz

## Future Work on Shape Determination

- Measurement data contain error
- better algorithm
- Choices of shape deviation parameters
- Extending the method to using frequencies, fields and external Qs where
- The forward problem is a complex nonlinear eigenvalue problem!
- Mesh smoothing (ITAPS)‏

- Meshes near
- pickup gap
- red: deformed
- black: original

## Linear and Nonlinear Eigensolvers

## RF Cavity Eigenvalue Problem 

- _Find frequency and field vector of normal modes:_

- “Maxwell’s Eqns in Frequency Domain”

## Cavity with Waveguide Coupling

- Vector wave equation with waveguide boundary conditions can be modeled by a **non-linear eigenvalue problem**

***Open***
***Cavity***

- Waveguide BC

- Waveguide BC

- Waveguide BC

*With *

- One waveguide mode per port only

**Notes:**
Closed cavity is an approximation for the real model. Usually the cavity has some sort of openings. For example, in accelearator cavity, we need put power into cavity and damp high-order-modes. These external coupling requires waveguide connected to cavity. Waveguide shape can be circular, rectangular or coax.

To model the open cavtity with waveguides, we can introduce a virtual BC on waveguide as shown. 

The problem is emerged from a cavity with external coupling through waveguides. We can put waveguide boundary conditions on the virtual boundaries. Curl-curl equation along with waveguide boundary conditions will be discretized using finite-element methods (Electric field E is expanded by a set of vector basis functions.)  and becomes a non-linear eigenvalue problem, where matrix K corresponds to curl-curl operator. M is mass matrix. Matrix W is waveguide matrix and has nonzero entries only on waveguide boundary.  Note that k_cj is a physical known quantity called cutoff.

## Cavity with Waveguide Coupling for Multiple Waveguide Modes

- Vector wave equation with waveguide boundary conditions can be modeled by a **non-linear eigenvalue problem (NEP)****‏**

*where*

**Notes:**
Closed cavity is an approximation for the real model. Usually the cavity has some sort of openings. For example, in accelearator cavity, we need put power into cavity and damp high-order-modes. These external coupling requires waveguide connected to cavity. Waveguide shape can be circular, rectangular or coax.

To model the open cavtity with waveguides, we can introduce a virtual BC on waveguide as shown. 

The problem is emerged from a cavity with external coupling through waveguides. We can put waveguide boundary conditions on the virtual boundaries. Curl-curl equation along with waveguide boundary conditions will be discretized using finite-element methods (Electric field E is expanded by a set of vector basis functions.)  and becomes a non-linear eigenvalue problem, where matrix K corresponds to curl-curl operator. M is mass matrix. Matrix W is waveguide matrix and has nonzero entries only on waveguide boundary.  Note that k_cj is a physical known quantity called cutoff.

## Physics Problems and Solver Options

**WSMP**

**MUMPS**

**SuperLU_Dist**

*Different solver options have different performance dynamics*

**Omega3P**

**Lossless**

**Lossy **
**Material**

**Periodic**
**Structure**

**External**
**Coupling**

**ESIL/with**
**Restart**

**ISIL w/ **
**refinement**

**Implicit****/Explicit**
** ****Restarted Arnoldi**

**SOAR**

**Self-Consistent**
**Iteration**

**Nonlinear**
**Arnoldi/JD**

**Physics Problems and ****Solver Options**

**Notes:**
Different physics, different eigensolvers, different linear solvers.”Complex” 

## Path to Simulate ILC RF Unit (3-cryomodule)‏

- Optimized ILC single cavity  routinely

- Simulated 4-cavity STF last year

- Simulating 8-cavity ILC Cryomodule this year

- Simulate ILC 3-cryomodule RF Unit
- - ~200M DOFs, further CS/AM advance needed, petascale

## Future Work for Eigensolvers	

- Parallelize AMLS, understand and improve its performance and scalability
- Nonlinear Jacobi-Davidson
- Choice of initial space
- Strategy for updating preconditioner and choice of preconditioners
- New algorithm development for NEP/LEP
- avoid shift-invert for interior eigenvalues
- LEP helps NEP (Self Consistent Iterations)

## Scalable Linear Solvers

## Linear Solver is Computational Kernel of Many Codes

- Indefinite Matrices
- Linear systems arising from shift-invert eigensolver in Omega3P
- Indefinite linear system from KKT conditions
- S-parameter computation in S3P
- Symmetric Positive Definite (SPD) Matrices
- From implicit time-stepping in T3P
- From thermal and mechanical analysis TEM3P
- From electro/magneto static analysis  Gun3P
- Issues in Petascale Electromagnetic simulations:
- Direct solver: memory usage, scalability of triangular solver
- Iterative solver: performance, effectiveness (preconditioner)

## Omega3P Scalability on Jaguar/XT with Iterative Linear Solver

- 1.5M tetrahedral elements   NDOFs = 9.6M  NNZ = 506M

- LCLS RF Gun

**Notes:**
Real problem instead of toy problem
Each rank has only about  < 400 Tetra when the problem ran on 4096 NCPUs!

## Scalability Using Sparse Direct Solver MUMPS

- Sparse Direct Solver is effective for highly indefinite matrices
- Scalability affected by performance of Triangular Solver

- Need more scalable Triangular Solvers

**Notes:**
SDS is very effective for solving linear system with highly indefinite matrices with multiple RHS
In Eigensolver, one needs to repeatedly solve linear systems with same matrices but different RHS. The Scalability of eigensolver is donimated by TS. 
This is performance of MUMPS TS for a matrix size of 2m. The n. of enetries in L is about 1 b. We can see that there is not much performance improvement as the ncpus increases.  On the other hand, the right hand side, is performance plot of another TS in PSPASES. We need a more scalable TS implementation.

## More “Memory-usage” Scalable Sparse Direct Solvers

- Maximal per-rank MU is 4-5 times than the average MU
- Once it cannot fit into Nprocs, it most likely will not fit into 2*Nprocs
- More “memory-usage” scalable solvers needed

- N=1.11M, nnz=46.1M Complex matrix

**Notes:**
People  pay more attention to speedup (performance scalablity) than memory-usage scalability. This is the per-rank memory usage of MUMPS for factorizeing 1.11 M complex matrix.  
Maximal per-rank MU 4-5 times
And it is not scalable. -- once 
Need a more “memory-usage” scalable” 

## Memory Saving Techniques 

- Single precision for factor matrix, iterative refinement to recover double precision accuracy (F)‏
- Domain-specific Preconditioners
- Factorize real part of the matrix (R)
- Real part is a good approximation to the complex matrix
- User single precision to factorize real part of the matrix (RF)
- Hierarchical preconditioners (FE order is the level) (HP)‏
- single precision for (1,1)-block (HPF)‏
- real part only for (1,1)-block (HPR)‏
- single precision & real part for (1,1)-block (HPRF)‏

## Testing Results for Complex Shifted Linear Systems

## Recent Progress of SuperLU(Xiaoye Li)

- Parallel symbolic factorization significantly reduces memory usage

- Matrix for DDS

- Matrix for ILC Cavity

## Future Work on Linear Solvers

- Direct versus iterative solvers, hybrid solvers
- Investigate applicability of out-of-core sparse direct solvers from TOPS
- Apply multigrid solvers from TOPS for SPD matrices
- Extend PSPASES to indefinite/complex matrices
- Develop more effective domain-specific preconditioners