Iterative methods with mixed-precision preconditioning for ill-conditioned linear systems in multiphase CFD simulations

A new mixed-precision preconditioner based on the iterative refinement (IR) method is developed for the preconditioned conjugate gradient (P-CG) solver and the multigrid preconditioned conjugate gradient (MGCG) solver in the multiphase thermal-hydraulic CFD code JUPITER [1]. In the IR preconditioner, linear systems are approximately solved using a block Jacobi method, which is optimized using fine block decomposition with 2D tiles to facilitate SIMD operations and software pipelining. Mixed-precision computing is designed to avoid influences of the roundoff errors in computing ill-conditioned matrices with extreme contrast of coefficients. Linear systems are normalized so that it can be computed within the dynamic range of FP16. All data is stored in FP16 to reduce memory access, while all computation is performed in FP32. The hybrid FP16/32 implementation keeps the similar convergence property as FP32, while the computational performance is close to FP16. The robustness of the FP16/32 implementation was demonstrated by scanning data formats with different bit lengths of the significand. The developed solvers are optimized on Fugaku (A64FX), and applied to ill-conditioned matrices with 90 billion DOFs in JUPITER. The P-CG and MGCG solvers with the new IR preconditioner show excellent strong scaling up to 8,000 CPUs, where 5.7 × and 3.4 × speedups are respectively achieved from the conventional solvers on Oakforest-PACS (KNL).