Advancing solver performance in large-scale computational fluid dynamics: generalizing the linelet preconditioner for the pressure correction equation

Advancing Solver Performance in Large-Scale Computational Fluid Dynamics: Generalizing the Linelet Preconditioner for the Pressure Correction Equation Author: Ramiro de Olazábal Supervisors: Dr. Oriol Lehmkuhl Barba Dr. Ricard Borrell A Thesis submitted to Universitat Politècnica de Catalunya for the degree of Doctor of Philosophy March 2025 Acknowledgements I want to express my sincere gratitude to my supervisors, Dr. Oriol Lehmkuhl and Dr. Ricard Borrell, for selecting me for this Ph.D. program. It was not an easy journey. Although the challenges, their guidance, support, and teachings have been invaluable, and I feel honored to have worked with them. I am deeply thankful to Dr. Herbert Owen for his support throughout this process. I am also thankful to Dr. Matias Avila for his patience with my many questions and Dr. Guillaume Houzeaux for smoothing my relationship with Alya. Special thanks go to Dr. Lucas Gasparino and Dr. Sarath Radhakrishnan for their support both within and outside BSC, to Carlos Arnedo for his patience while teaching me to trace, and to Samuel Gomez and Marc Silvestre i Claros for helping me with ANSA every time I needed it. I sincerely thank all my colleagues and friends at BSC; this journey would have been much more challenging without them. I would like to acknowledge the Barcelona Supercomputing Center for providing the resources and facilities necessary for conducting this research. I am also grateful for the financial support from Spain’s Ministerio de Ciencia e Innovación through the grant “Ayudas para contratos predoctorales para la formación de doctores” (Ref: PRE2018-086548), which allowed me to advance my studies and research. There is no way to thank my family enough. I am especially grateful to my parents and sister, who have always been by my side despite the distance, offering unwavering support and encouragement. I must also express my gratitude to Dr. Fernando Lombardo. I owe him the opportunity to work in such an incredible and intellectually challenging place as BSC. Your support and belief in me have been genuinely decisive in my career. I am equally thankful to the Echavarría and Belmonte families for warmly welcoming me when I arrived in Barcelona and making me feel at home. Their kindness has greatly eased the challenges of these years. I cannot forget my friends—both lifelong companions and those I made over these years. Thank you all: to David Oks, for easing my transition into Bari celona; to Santi, for the engaging conversations and shared beers; to Cristóbal, for the stimulating talks and wines; to Pablo, Bruno, and David, for the trips and good times; to Nico, for the shared moments; and Emi and Juani, for bringing a touch of Argentina’s warmth. I am deeply grateful to everyone who contributed to this thesis in one way or another for their support, guidance, and encouragement. ii Abstract The efficient solution of the Poisson equation is a fundamental challenge in Computational Fluid Dynamics (CFD), particularly in large-scale simulations involving high-Reynolds-number flows. The strong anisotropy introduced by thin boundary layers requires highly stretched elements, which, in turn, significantly impact the convergence of iterative solvers. In such cases, the Preconditioned Conjugate Gradient (PCG) solver, commonly used for pressure correction, exhibits slow convergence due to the deteriorated conditioning of the system matrix. To address this issue, the Linelet Preconditioner has been widely adopted, as it exploits the underlying mesh anisotropy by constructing linelets—one-dimensional structures aligned with the strongest couplings—and applying specialized operations along these segments. However, traditional implementations of the Linelet Preconditioner operate within each domain partition independently, resulting in degraded performance as the number of partitions increases in massively parallel simulations. This thesis presents the Global Linelet Preconditioner (GLP), a preconditioning strategy designed to overcome these limitations and improve scalability in extreme-scale CFD applications. The method extends and generalizes the traditional linelet approach by introducing a communication step within the preconditioning operation, allowing interdomain coupling and preserving connectivity across partition boundaries. This modification significantly enhances convergence rates in highly anisotropic meshes by ensuring that the strongest couplings in the linear system are treated effectively, regardless of domain decomposition constraints. A key contribution of this work is the development of a purely algebraic linelet construction algorithm, which eliminates the need for geometric information when defining linelets. While conventional methods rely on explicit mesh structures to determine anisotropic directions, the algebraic approach constructs linelets based solely on matrix properties, allowing greater flexibility and applicability to general unstructured meshes. Furthermore, the geometricbased construction was also explored and integrated within the framework, demonstrating superior performance in structured meshes with well-defined anisotropic features. The comparison between the geometric and algebraic approaches revealed that while the former achieves better performance when clear directional stiffness is present, the latter provides a robust alternative when mesh topology is complex or unavailable. iii The effectiveness of GLP was assessed through extensive numerical experiments, including benchmark problems and real-world CFD applications such as the 30P30N high-lift airfoil, the Stanford diffuser, and the DrivAer model. Results demonstrated that GLP significantly improves solver convergence over existing preconditioners, including previous versions of the linelet preconditioner, particularly in cases where a high percentage of elements lie within the boundary layer. Performance analyses revealed that while GLP incurs a higher preprocessing cost due to linelet construction and communication, these overheads are outweighed by the substantial reduction in PCG iterations, leading to overall computational savings in large-scale simulations. In addition to improving convergence, GLP introduces a partition-agnostic formulation, making it independent of the domain decomposition strategy. Unlike traditional preconditioners, which are sensitive to mesh partitioning, GLP maintains its numerical performance across varying decomposition configurations, enabling more flexible and balanced parallel execution. The parallel implementation of the method, tailored for High Performance Computing (HPC) environments, ensures scalability across a wide range of core counts, as demonstrated by detailed scalability analyses. Keywords: Computational Fluid Dynamics, Poisson Equation, Preconditioned Conjugate Gradient, Linelet Preconditioner, High-Performance Computing, Anisotropic Meshes, Parallel Solvers, Algebraic Preconditioners, Domain Decomposition. iv Contents Acknowledgements i Abstract iii Contents v List of Figures vii List of Tables xii List of acronyms and abbreviations xiv 1 Introduction and fundamentals 1 1.1 Motivation............................. 1 1.2 Objectives and thesis structure . . . . . . . . . . . . . . . . . . 4 1.3 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . 5 1.3.1 Reference systems . . . . . . . . . . . . . . . . . . . . . 6 1.3.2 Fluid mechanics . . . . . . . . . . . . . . . . . . . . . . 8 Complete Set of Equations . . . . . . . . . . . . . . . . 8 1.4 Discretization of the Navier-Stokes equations . . . . . . . . . . 9 1.4.1 Space Discretization . . . . . . . . . . . . . . . . . . . 10 Weak Formulation . . . . . . . . . . . . . . . . . . . . 10 Weak Formulation of the Navier-Stokes Equations . . . 11 Time-dependent formulation: . . . . . . . . . . . 11 Finite Element Method . . . . . . . . . . . . . . . . . . 12 Finite Element Discretization of the Navier-Stokes Equations.................... 14 1.4.2 Time Discretization . . . . . . . . . . . . . . . . . . . . 15 Fractional Step Schemes for the Navier-Stokes Equations 16 1.5 Methods for Solving Linear Systems . . . . . . . . . . . . . . . 17 1.5.1 Direct Methods . . . . . . . . . . . . . . . . . . . . . . 18 1.5.2 Multigrid Methods . . . . . . . . . . . . . . . . . . . . 18 1.5.3 Iterative Methods . . . . . . . . . . . . . . . . . . . . . 18 KrylovMethods...................... 19 Conjugate Gradient Method . . . . . . . . . . . . . . . 19 Preconditioned Conjugate Gradient Method . . . . . . 20 1.5.4 Different types of preconditioners . . . . . . . . . . . . 22 v Contents 1.6 Linelet preconditioner . . . . . . . . . . . . . . . . . . . . . . . 23 1.6.1 Caracterization of the linelet preconditioner . . . . . . 25 1.6.2 Effect of the integration rule on the linelet preconditioner 27 1.6.3 Parallelization of the linelet preconditioner . . . . . . . 29 1.7 Auxiliary methods . . . . . . . . . . . . . . . . . . . . . . . . . 33 1.7.1 TDMA: Tridiagonal Matrix Algorithm . . . . . . . . . 33 1.7.2 The Schur-complement method . . . . . . . . . . . . . 33 1.8 Conclusions ............................ 35 References................................ 35 2 Efficient Parallelization of the Tridiagonal Matrix Algorithm 39 2.1 Various approaches . . . . . . . . . . . . . . . . . . . . . . . . 40 2.1.1 Previous approaches . . . . . . . . . . . . . . . . . . . 40 2.1.2 Our approach: agnostic to the domain partition . . . . 40 2.2 Parallelization of the TDMA using the Schur-Complement Method............................... 41 2.2.1 An efficient way to compute each slice’s contribution tother.h.s. ........................ 41 Direct computation of the product As,pA−1 p,p ...... 42 2.2.2 Parallel approach to solve the linear system . . . . . . 45 2.3 Conclusions ............................ 49 References................................ 49 3 Assembling the preconditioning matrix 51 3.1 Assembling the slices . . . . . . . . . . . . . . . . . . . . . . . 52 3.1.1 Algebraic slice assembly within each subdomain . . . . 52 Remarks on choosing the next node in the slice . . . . 53 3.1.2 Geometric slice assembly within each subdomain . . . 58 3.2 Halocleaning ........................... 59 3.2.1 Each process has the global ID of its nodes and the nodes to which they are coupled via the halo . . . . . . 61 3.2.2 Each process does not have the global ID of its nodes and the nodes to which they are coupled via the halo . 62 3.3 Growing the linelets: joining slices from different subdomains 66 3.3.1 Selecting a specific set of slices to be the starting point of their respective linelets . . . . . . . . . . . . . . . . 67 3.3.2 Growing the linelets without specifying a starting point 70 Loop finding strategies . . . . . . . . . . . . . . . . . . 76 Loop cleaning strategy . . . . . . . . . . . . . . . . . . 78 Growing the linelets . . . . . . . . . . . . . . . . . . . . 79 3.4 Conclusions ............................ 79 References................................ 80 4 Implementation and performance 81 4.1 Implementation .......................... 82 4.2 Workflow.............................. 82 4.3 Computational Resources . . . . . . . . . . . . . . . . . . . . . 82 4.4 Numerical verification of the solver . . . . . . . . . . . . . . . 83 4.5 Analysis of performance . . . . . . . . . . . . . . . . . . . . . 98 4.5.1 Execution Analysis . . . . . . . . . . . . . . . . . . . . 98 vi 4.5.2 Scalability Analysis and Performance Metrics . . . . . 101 4.6 Preprocessing stage . . . . . . . . . . . . . . . . . . . . . . . . 106 4.7 Conclusions ............................ 109 References................................ 113 5 Numerical Applications – Solving Real-World Scenarios 114 5.1 30P30NAirfoil........................... 115 5.1.1 Computation and Communication Time Analysis . . . 116 5.1.2 Computational Cost per PCG Iteration . . . . . . . . . 117 5.2 Stanford Diffuser . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.3 DrivAerModel........................... 127 5.4 Performance Analysis . . . . . . . . . . . . . . . . . . . . . . . 130 5.4.1 Preprocessing Stage . . . . . . . . . . . . . . . . . . . . 130 5.4.2 Preconditioning Step . . . . . . . . . . . . . . . . . . . 132 5.5 Scalability Analysis and Performance Metrics . . . . . . . . . 134 5.5.1 SpeedUp ......................... 134 5.5.2 Communication Efficiency . . . . . . . . . . . . . . . . 135 5.5.3 Load Balance . . . . . . . . . . . . . . . . . . . . . . . 136 5.5.4 Parallel Efficiency . . . . . . . . . . . . . . . . . . . . . 137 5.6 Conclusions ............................ 139 References................................ 141 6 Conclusion 142 6.1 Limitations............................. 143 6.2 FutureWork............................ 145 References................................ 145 Appendices 146 A List of publications and contributions 147 A.1 Academic writings . . . . . . . . . . . . . . . . . . . . . . . . . 147 A.2 Conferences and workshops . . . . . . . . . . . . . . . . . . . . 147 B Detailed description on growing the linelets 148 B.1 Growing the linelets from selected slices . . . . . . . . . . . . 148 B.2 Growing the linelets without specifying a starting point . . . . 153 List of Figures 1.1 Eulerian, Lagrangian, and ALE reference systems. . . . . . . . . . 7 vii List of Figures 1.2 Linelet structure derived from the strongest couplings in the mesh and the tridiagonal matrix obtained by exclusively considering these dominant connections. . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.3 Comparison of the number of iterations required to reach a tolerance of 10 −6 for the linelet and diagonal preconditioner in the PCG Method on meshes with varying maximum aspect ratio and boundary layer coverage. . . . . . . . . . . . . . . . . . . . . . 26 1.4 Convergence of the residual for the linelet preconditioner for different integration rules. . . . . . . . . . . . . . . . . . . . . . . . 28 1.5 Convergence of the preconditioned residual for the linelet preconditioner for different integration rules. . . . . . . . . . . . . . . . . 28 1.6 Same mesh as depicted in Figure 1.2a, but distributed across five distinct subdomains. The conventional method involves restricting linelets to a single domain partition, excluding inter-subdomain couplings from the preconditioning matrix. . . . . . . . . . . . . . 30 1.8 Convergence of the PCG method for the case of a linelet preconditioner applied to a partitioned mesh. . . . . . . . . . . . . 30 1.7 Example of a mesh for the numerical simulation of an airplane wing. On the right, a closer view provides detailed insight into the prismatic boundary layer. Typically, a significant percentage of the mesh elements are located here. The pronounced anisotropy in this region introduces stiffness to the discrete Poisson’s equation. . . . 31 1.9 Number of iterations required for the case shown in Figure 1.7 to reach a tolerance of 10−3at each time step. . . . . . . . . . . . . . 32 1.10 Number of iterations required for the case shown in Figure 1.8 to reach a tolerance of 10−6at each time step. . . . . . . . . . . . . . 32 2.1 Definition of the terms "slice" and "linelet": We introduce the term "slice" to refer to the local structure of nodes that are connected within a single subdomain. In contrast, "linelet" refers to the global structure arising from many slices’ connection. Interface nodes are depicted as empty nodes, and inner nodes as filled ones. . . . . . . 41 2.2 Diagram showing the tridiagonal structure of a linelet preconditioningmatrix. ............................. 43 2.3 Structure of a reordered linelet preconditioning matrix according to the Schur algorithm. . . . . . . . . . . . . . . . . . . . . . . . . 44 2.4 Schematic representation of the proposed parallel algorithm to solve the tridiagonal preconditioning matrix. Each linelet is assigned to a specific process responsible for solving the linear system at the interface. ................................ 48 3.1 Iterative process to build slices. . . . . . . . . . . . . . . . . . . . 54 3.2 Sketch showing nodes and couplings for an example mesh. Nodes h,iand jshould belong to the same slice. . . . . . . . . . . . . . . 55 viii CHAPTER 1 Introduction and fundamentals 1.1 Motivation Many physical laws are expressed through differential equations, as they provide a fundamental framework for describing the dynamics of physical systems. However, analytical solutions exist only for a limited set of simplified cases, necessitating numerical methods for more realistic scenarios. For example, in Computational Fluid Dynamics (CFD), accurate simulations require solving large linear systems derived from the Navier-Stokes equations, which govern fluid motion through nonlinear, coupled partial differential equations. Despite significant progress in numerical methods, solving these equations remains challenging. Moreover, proving the existence and smoothness of solutions for the Navier-Stokes equations is recognized as one of the millennium problems in mathematics [5]. Numerical solutions provide insights into critical aspects of fluid behavior, such as turbulence, heat transfer, and pressure distribution, which are essential for optimizing designs, ensuring safety, and advancing technology. Consequently, CFD has become an indispensable tool across various industries, including aerospace, automotive, energy, medicine, and environmental science. Applications range from designing more aerodynamically efficient aircraft and reducing vehicle emissions to optimizing wind turbines and modeling blood flow in medical research. While numerical methods for solving the Navier-Stokes system of equations governing fluid flow were proposed as early as the 19th century, the advent of digital computing in the mid-20th century enabled more complex algorithms to be developed and applied to more realistic scenarios. One of the main challenges in finding an accurate solution to this set of equations is the rise of turbulent behavior in the flow, where the velocity field becomes chaotic and unpredictable. This phenomenon is characterized by a wide range of interacting scales in which the energy cascades from large, unstable eddies to smaller ones until it is dissipated into heat. The fact that the NavierStokes equations give rise to this chaotic and multiscale behavior makes the direct solution of these equations infeasible for most turbulent flows. The growth of massively parallel computing architectures in recent decades has significantly expanded the applicability of these methods to large-scale problems. 1 1.1. Motivation For incompressible flows, the fractional step projection method is commonly used to decouple the velocity and pressure fields while ensuring mass conservation [22]. This approach requires solving a Poisson equation for pressure correction at each time step, which is often the most computationally expensive component of the solver and also one of the most difficult to parallelize. In our experience, we have observed that this equation must be solved with the highest possible precision. The reason behind this stringent accuracy criterion lies in the necessity to uphold mass conservation. When solved with low precision thresholds, the conservation of mass tends to be compromised, leading to non-convergence of the problem. Also, the time taken to solve this Poisson equation grows with the time step employed, making efficient solvers essential for large-scale simulations. Our observations show that it consumes between 40% and 60% of the solver time. Moreover, for highly anisotropic meshes, such anisotropy degrades the conditioning of the linear system associated with the Poisson equation. Numerical methods for solving large linear systems can be classified into direct and iterative solvers. While direct methods provide exact solutions, they are impractical for large-scale problems due to their O ( n2 )computational and memory complexity. Iterative methods, such as Krylov subspace solvers and multigrid methods, are more efficient for large systems. However, their convergence rate is highly dependent on the conditioning of the system matrix and the choice of preconditioners. Preconditioners play a crucial role in accelerating the convergence of iterative solvers by improving the numerical properties of the system matrix. In CFD, the presence of highly anisotropic meshes—particularly in boundary layers—leads to poorly conditioned linear systems, which slow down convergence. Specialized preconditioners are required to mitigate these issues and ensure efficient simulations. A critical challenge in CFD is modeling turbulence in the flow. Turbulence is characterized by chaotic fluid motion in which a wide range of spatial and temporal scales interact by mixing, heat transfer, and momentum transport. In high-Reynolds-number flows, turbulent fluctuations span from large energy-containing eddies to small dissipative scales, making their full resolution computationally prohibitive. As the energy cascades down through a hierarchy of scales, it is eventually dissipated by viscous forces at the smallest scales, —known as the Kolmogorov scales. This phenomenon is called the Kolmogorov cascade, which describes how energy is transferred from large to progressively smaller scales. The complexity of this phenomenon has motivated the development of various modeling strategies that aim to capture the essential features of turbulence without resolving every scale. In practical CFD simulations, accurately capturing the complex nature of turbulence requires adopting modeling strategies, each with pros and cons regarding fidelity and computational cost. Direct Numerical Simulation (DNS) resolves all turbulent scales but is limited to low Reynolds numbers due to its high computational cost. Large Eddy Simulation (LES), on the other hand, resolves the large-scale eddies while modeling the effects of smaller scales, 2 1.1. Motivation balancing fidelity and efficiency. Reynolds-Averaged Navier-Stokes (RANS) methods further reduce the computational cost by averaging the turbulent fluctuations. Additionally, LES often employs wall modeling techniques to accurately capture near-wall behavior without requiring prohibitively fine meshes in these regions. These modeling choices directly affect the numerical properties of the discretized equations. For instance, employing Wall-modelled Large Eddy Simulation (WMLES) in high-Reynolds-number flows can lead to anisotropic meshes within the boundary layer, which can exacerbate the ill-conditioning of the linear systems. This scenario underscores the need for robust preconditioning strategies capable of handling the stiffness introduced by steep gradients and stretched elements. A related and equally critical challenge in fluid simulations is the accurate resolution of the boundary layer, a thin region near solid surfaces where velocity gradients are steep. The physical characteristics of the boundary layer are influenced by the Reynolds number, which indicates whether the flow is laminar (smooth) or turbulent (chaotic). In this region, viscosity effects are highly significant, and the resolution of this region requires highly anisotropic grids. The flow in this region remains predominantly laminar for low Reynolds numbers. In contrast, it transitions from laminar to turbulent at a certain distance from the leading edge at high Reynolds numbers. Although laminar boundary layers generate less friction, they are unstable and can quickly become turbulent. Properly simulating the transition from laminar to turbulent flow is crucial for accurately predicting aerodynamic performance and minimizing energy losses caused by flow separation. In high Reynolds number flows, boundary layers are subject to significant anisotropy due to steep velocity gradients near surfaces, necessitating fine grid resolution to capture these effects accurately. This results in highly stretched elements in the mesh, which contribute to the stiffness and ill-conditioning of the system matrix. For such problems, standard preconditioners like diagonal, ILU, and multigrid methods may not perform optimally. It has been observed that while multigrid methods can reduce the number of PCG iterations, their total computational cost remains comparable to that of simpler preconditioners [30]. Furthermore, they struggle with highly stretched elements, making them less effective for boundary-layer-dominated problems. To address the challenges posed by anisotropic grids, the linelet preconditioner was introduced as a specialized approach that approximates the strongest couplings in the system through a set of decoupled 1D problems. Previous studies [23,28,1,4,12,18,20,9,21] have demonstrated its effectiveness, showing significant improvements over conventional preconditioners such as diagonal, ILU, and multigrid [23,28]. In particular, [28] reported that the linelet preconditioner reduces the number of PCG iterations by at least 50% and decreases total CPU time by up to 40% compared to diagonal [24], multigrid 3 1.2. Objectives and thesis structure [32], and ILU [24] preconditioners. However, traditional implementations impose a critical constraint: linelets must be confined within a single subdomain, limiting scalability in parallel computations. Existing approaches attempt to mitigate this by modifying either the linelet structure or the partitioning algorithm to ensure that each linelet remains contained within a subdomain [23]. This restriction, however, reduces flexibility and hinders the preconditioner’s efficiency in large-scale simulations. Despite the widespread use of linelet-based preconditioners, no existing solver allows linelets to span multiple processes in a distributed computing environment. This work addresses these limitations by developing a purely algebraic linelet preconditioner that operates directly on the system matrix, independent of mesh geometry. This new approach enables linelets to extend across multiple processes, improving scalability and performance in parallel computing environments. Unlike geometric-based methods, this purely algebraic formulation preserves the key advantages of linelet preconditioning while eliminating the constraints imposed by domain partitioning. Finally, the increasing complexity of numerical simulations in CFD aligns with broader trends in High-Performance Computing (HPC). While traditional performance improvements have followed Moore’s law, future advancements will depend on innovations at higher levels of the computing stack, including hardware architecture, software optimization, and algorithmic development. This transition requires efficient numerical solvers that can fully exploit modern computing capabilities. Developing scalable, high-performance preconditioners is essential for advancing large-scale CFD simulations. By improving solver efficiency, this work contributes to reducing computational costs and enabling more detailed and accurate simulations, ultimately supporting advancements in aerospace engineering, energy systems, environmental modeling, and other critical fields. 1.2 Objectives and thesis structure The primary objective of this thesis is to develop an efficient preconditioning strategy for large-scale CFD simulations, with a particular focus on overcoming the challenges posed by highly anisotropic meshes. Simulating high-Reynoldsnumber turbulent flows is frequently hindered by the slow convergence of iterative Poisson solvers, mainly due to the ill-conditioning caused by thin boundary layers and severe mesh stretching. Building on the work of Soto, Löhner, and Camelli [28], we extend the concept of the linelet preconditioner into a global formulation. This Global Linelet Preconditioner (GLP) is designed to preserve the interdomain couplings in distributed-memory environments, thereby accelerating convergence while remaining robust in massively parallel simulations. Our approach integrates both algebraic and geometric perspectives: the algebraic formulation allows for application without direct reliance on 4 1.3. Governing Equations geometric information, while the geometric construction takes into account directional features of the mesh to better capture anisotropy. In essence, the core objective of this thesis can be translated into two deceptively simple yet challenging questions. First, how can one construct a tridiagonal matrix that is an effective preconditioner for the Poisson equation in a distributed-memory environment while remaining agnostic to the domain partition? Second, how can we efficiently parallelize the Tridiagonal Matrix Algorithm (TDMA) algorithm when the matrix entries are spread over many processes? The answers to these questions lie at the core of our work. In particular, the solution to the first question generalizes the construction approach presented in [28] by incorporating interdomain couplings in addition to the local ones. The thesis is organized to introduce the reader to the theoretical foundations and then to the practical implementation and validation of our method. Chapter 1 establishes the fundamental principles of fluid mechanics and the numerical treatment of the Navier–Stokes equations, emphasizing the computational challenges associated with solving these equations on highly anisotropic meshes. Chapter 2 addresses the efficient parallelization of the TDMA by incorporating the Schur complement method, the chosen approach for maintaining interdomain coupling across different processes. In Chapter 3, we detail the assembly process of the preconditioning matrix: local linelet slices are constructed within each subdomain and then joined across subdomain boundaries to form a coherent global structure. Chapter 4 focuses on the implementation within the Alya simulation software and presents an extensive performance evaluation, including analyses of scalability, load balance, and communication efficiency. Finally, Chapter 5 applies the GLP preconditioner to real-world CFD scenarios—the 30P30N high-lift airfoil, the Stanford diffuser, and the DrivAer model—demonstrating its effectiveness in practical, large-scale simulations. Finally, the conclusion not only summarizes the achievements and key contributions of the work but also discusses its limitations and outlines directions for future research. 1.3 Governing Equations This section presents some of the equations and numerical methods usually employed in CFD. First, we introduce the Navier-Stokes equations, which describe fluid motion. Then, the finite element method is introduced for spatial discretization, providing the framework for approximating the governing equations on complex geometries. Then, the fractional step projection method is used to propagate the solution in time. This method decouples the pressure and velocity calculations while maintaining the incompressibility condition. This method involves solving a Poisson equation for pressure correction at each time step. The resolution of this Poisson equation plays a critical role in maintaining mass conservation and ensuring convergence. Solving the pressure Poisson equation becomes increasingly challenging on 5 1.3. Governing Equations highly anisotropic meshes, as these can lead to stiffness and poor conditioning in the system matrix. To address this issue, preconditioners are used to improve the convergence rate of iterative solvers by reducing the condition number of the matrix, making the iterative solvers more efficient when handling large linear systems. This ensures computational efficiency and enables the effective resolution of such systems, which is a central focus of this study. 1.3.1 Reference systems Let us introduce the theoretical framework by briefly discussing the reference systems commonly used in fluid mechanics. Indeed, before attempting to describe the physics of a given fluid, it is necessary to define the reference system—or point of view—that will be used to describe its dynamics. In this regard, three types of reference systems are commonly used in fluid mechanics. First, we have the Lagrangian reference system (also known as the current or material description), which we will denote as RL . In this system, the observer follows a specific fluid particle as it moves through time and space. Thus, the position and velocity of a particle that at a given time t0 is located at position x0will be given by x=x(t, x0) U(t, x0) = ∂x(t, x0) ∂t (1.1) The rest of the relevant quantities are described in the same way, such as the volume differential δV ( t0,x0 )of the fluid element that at time t0 is in the surroundings of x0 , its mass δM ( t0,x0 ), density ρ ( t0,x0 ) = δM ( t0,x0 ) /δV ( t0,x0 ), temperature T ( t0,x0 ), etc. Thus, from this perspective, describing a particular property of a continuous medium involves expressing it as a function of all initial positions x0and all required time instants t:A(t0,x0). On the other hand, the Eulerian description (also known as spatial configuration) is carried out in a reference system fixed to the laboratory, which we will denote as RE . Unlike the Lagrangian description, where a fluid particle is followed in its motion and its fluid properties are described in terms of the property of the particle at time t0 located at position x0 , in the Eulerian description, we describe the property of the fluid particle that is located at position x at time t . This will, therefore, be a field description of the fluid. Thus, the velocity field will be given by U ( t, x ), the density by ρ ( t, x ), the temperature by T ( t, x ), and similarly, for any other property defined on the fluid a ( t, x ). Finally, there is a third option, ALE (referential configuration), which we will denote as RA . In this system, the fluid properties are described in an independent reference system, typically fixed to the computational grid. This is particularly relevant when working with mesh distortions, as it allows for preserving an accurate description of interfaces under significant mesh distortions. 6 1.3. Governing Equations Since the three representations are different ways of visualizing the same physical phenomenon, it is possible to transition from one to another through the appropriate transformations. Thus, we can consider the family of invertible mappings defined by φ:RE×[t0, tend)→RA×[t0, tend) Ψ : RA×[t0, tend)→RL×[t0, tend) Φ : RL×[t0, tend)→RE×[t0, tend) (1.2) Figure 1.1: Eulerian, Lagrangian, and ALE reference systems and their transformations. Since the three descriptions are equivalent, if, for example, the Lagrangian description of a fluid property is known, the Eulerian description will be given by: a(t, x) = A(t, x0(t, x)),(1.3) where x0(t, x)is the equation that results from solving for x0in the equation x(t, x0) = x,(1.4) which relates the original position x0 that a given particle had at time t0 with the position it occupies at time t . Due to particle identity, the previous relation will always be bijective and, therefore, invertible. Similarly, given the Eulerian description of a medium, the corresponding Lagrangian description will be given by A(t, x0) = a(t, x(t, x0)),(1.5) where the following equation must be solved beforehand x(t, x0) = x0+Zt t0 u(τ, x(τ, x0))dτ, (1.6) which results from integrating the equation 7 1.3. Governing Equations ∂x(t, x0) ∂t =U(t, x0) = u(t, x0).(1.7) Given that we will not work with mesh movement in this thesis, the ALE description will coincide with the Eulerian one. We will use the latter throughout this work. 1.3.2 Fluid mechanics Since this work aims to improve the efficiency of solvers for the Poisson equation, we consider it important to present a detailed background and the role the Poisson equation plays in CFD. For that, we present here the complete set of Navier-Stokes equations. This will allow us to understand not only the discretization methods used but also the role that the Poisson equation plays in their context and the importance of having fast and efficient methods for its solution. We will also consider the case of a Newtonian incompressible fluid. This approximation, as we will see, is valid for low Mach number flows, that is, flows in which the velocity of the fluid is much smaller than the speed of sound. In this section, we formulate the Navier-Stokes equations in the Eulerian frame of reference. Complete Set of Equations The Navier-Stokes equations describe the motion of fluid substances under the influence of various forces. For an incompressible fluid moving in a domain Ω bounded by Γ = ∂ Ωover the time interval ( t0, tf ), the governing equations are given by ρ∂tui+ρuj∂jui+∂ip−∂j(µ(∂jui+∂iuj)) = fion Ω×(t0, tf),(1.8a) ∂juj= 0 on Ω×(t0, tf).(1.8b) for i, j = 1 , ..., d , with d being the dimension of the domain. Here, ui is the velocity field, p is the pressure, ρ is the fluid density, µ is the dynamic viscosity, and fi is the momentum source. The first equation in (1.8) represents the conservation of momentum, while the second equation enforces incompressibility. To fully describe the fluid motion, the system in (1.8) must be supplemented with appropriate boundary conditions over Ω × ( t0, tf ). Defining Γ = ∂ Ω, we distinguish between Dirichlet conditions, which prescribe velocity, and Neumann conditions, which prescribe traction. Let Γ D and Γ N be the portions of the boundary where these conditions apply, satisfying Γ D∪ Γ N = Γ and Γ D∩ Γ N = ∅ . The boundary conditions are given by ui=uD ion ΓD×(t0, tf),(1.9a) σijnj=tN ion ΓN×(t0, tf),(1.9b) where nj is the unit normal to the boundary Γ N , and σij is the Cauchy stress tensor: 8 1.4. Discretization of the Navier-Stokes equations σij =−pδij +µ(∂iuj+∂jui).(1.10) Here, tN irepresents the imposed traction at the Neumann boundaries. Since the equations in (1.8) contain first-order temporal derivatives of velocity, an initial condition must also be imposed: ui(x, t0) = u0 i(x)on Ω.(1.11) If the fluid density ρ and viscosity coefficient µ are assumed constant, introducing the kinematic viscosity ν = µ/ρ and redefining p to be the pressure per unit density (p→p/ρ), the system simplifies to ∂tui+uj∂jui+∂ip−ν∂j(∂jui+∂iuj) = fion Ω×(t0, tf),(1.12a) ∂juj= 0 on Ω×(t0, tf),(1.12b) ui=uD ion ΓD×(t0, tf),(1.12c) σijnj=tN ion ΓN×(t0, tf),(1.12d) ui(x, t0) = u0 i(x)on Ω× {t0}.(1.12e) for i, j = 1, ..., d. 1.4 Discretization of the Navier-Stokes equations This section focuses on the discretization of the Navier-Stokes equations, which converts them from a continuous mathematical model into a form that can be solved using computational methods. We begin by exploring various spatial discretization techniques suited to different problem characteristics. Among the most common methods, we chose the FEM, which is well-suited for complex geometries. Next, we introduce the weak formulation of the Navier-Stokes equations, a less restrictive approach than the strong form. Thus enabling solutions in broader function spaces. This foundation allows the application of the FEM to solve the equations, including time-dependent problems, by incorporating suitable function spaces for temporal solutions. In FEM, the solution is approximated as a linear combination of shape functions defined on the mesh nodes, reducing the problem of determining the nodal unknowns. We then move into the time discretization process for the Navier-Stokes equations. The discretized form of the equations, incorporating both velocity and pressure, is introduced next. The goal is to find Un+1 and Pn+1 , which are the nodal values of velocity and pressure at the subsequent time step. The system of equations is outlined, setting the stage for the fractional step schemes. At this point, fractional step schemes are presented as an approach to decouple velocity and pressure in the time discretization. These schemes simplify the problem into separate steps, making the solution process more efficient. The system is solved in stages, with the momentum advanced explicitly and a Poisson equation solved for pressure stabilization. The fractional step method leads to a Poisson-like equation for pressure, which requires the solution of a linear system of the form 9 1.4. Discretization of the Navier-Stokes equations L ∆ P = b , where L approximates the Laplacian operator. This step is one of the most time-consuming and challenging to parallelize, and thus, achieving a solution is expensive in time and computational resources: solving it efficiently is crucial for the performance of the numerical method. Finally, we highlight how element dimensions affect the overall matrix Lcontributions. 1.4.1 Space Discretization The choice of discretization method depends on the complexity of the problem, the geometry of the domain, and the desired properties of the solution. The Finite Elements Method (FEM) constructs local approximations of the solution using local data, assembling them into a global approximation. One of its main advantages is its ability to handle complex, irregular geometries due to its use of flexible, unstructured meshes that conform well to intricate boundaries. Unlike the Finite Difference Method (FDM), which relies on structured grids and is less suited for irregular domains, FEM allows for local refinement and adaptivity, making it particularly useful in regions with steep gradients or high solution variability. On the other hand, the Finite Volumes Method (FVM) evaluates exact expressions for the average value of the solution over a control volume. It defines the solution as cell-averaged values, requiring reconstruction or interpolation methods to approximate values at boundaries, typically using neighboring cell data. Unlike FEM, which supports higher-order basis functions without the need for reconstruction, FVM enforces local flux conservation, ensuring that the net flux across each control volume is balanced. In contrast, FEM does not guarantee local conservation but ensures global conservation of fluxes. Given its widespread adoption in industry and its suitability for handling multiphysics problems, this work focuses on the FEM for spatial discretization. However, it is important to note that the methods developed in this thesis can be extended to other discretization approaches. Weak Formulation Before applying the FEM to the Navier-Stokes equations, it is necessary to introduce the weak formulation of a PDE. This formulation is the basis of the Finite Element Method. The main advantage of this approach is that it is less restrictive and relaxes some of the original PDE constraints. Consider a differential equation for a function u(x)in a domain Ω: L {u(x)}=f(x),(1.13) with the appropriate boundary conditions. Equation (1.13) is known as the strong form of the PDE. We define the function space V = H1 (Ω) as the space of functions v such that v∈L2 (Ω) and ∂iv∈L2 (Ω), where L2 (Ω) is the space of square-integrable functions, and H1 (Ω) consists of functions whose first derivatives also belong to L2 (Ω). Multiplying equation (1.13) by a test function v∈V and integrating 10 1.5. Methods for Solving Linear Systems the projection step consists of solving a Poisson equation for the pressure stabilization [7]. This means that at this step, a linear system of the form L∆P=b, (1.36) needs to be solved, where matrix L is the approximation to the Laplacian operator, and ∆ P is the increment in pressure (which is the unknown to be solved for). Note also that only the right-hand side (r.h.s.) b is affected by the properties of the flow, whereas L depends only on the properties of the mesh. This associated algebraic system is the implicit part of the time step and is one of the most time-consuming and difficult to parallelize. Therefore, choosing a suitable numerical method is crucial to obtaining a solution without being expensive in terms of computational resources. This matrix plays a fundamental role in many applications, such as heat transfer, structural mechanics, and computational fluid dynamics. Hence, it is important to have an efficient method to compute and solve it. In the following sections of this work, we will express equation (1.36) in the form Ax = b , where A represents the matrix that approximates the Laplacian operator and x denotes the pressure increment. This notation helps streamline our presentation and highlights that the proposed method is applicable for solving equation (1.36) and any linear system that follows the linelet structure, which will be described in the subsequent sections. 1.5 Methods for Solving Linear Systems The solution of large linear systems can be addressed using either direct or iterative techniques. Direct methods, such as Gaussian elimination and factorization (e.g., LU or Cholesky), theoretically provide exact solutions after a finite number of operations. However, they are impractical for large-scale problems due to their high computational and memory demands. Specifically, direct solvers require O ( n3 )arithmetic operations and O ( n2 )memory storage, where n is the number of unknowns [11,17]. These requirements make direct methods unsuitable for large-scale applications. Additionally, floating-point errors can accumulate, preventing exact solutions in practice. On the other hand, Iterative methods are particularly advantageous for sparse systems, where the scalability issues of direct methods become prohibitive as the problem size n grows. Iterative solvers generally require only O ( n ) arithmetic operations and O ( n )memory, making them preferable for large problems. These methods start from an initial guess and refine it iteratively until the residual falls below a predefined threshold, ensuring sufficient accuracy while avoiding unnecessary computations. The convergence of iterative methods depends on the conditioning of the matrix and the quality of the initial guess, and they are widely used in practice due to their adaptability and efficiency. Among iterative techniques, three main categories exist: stationary methods, multigrid methods [3,31], and Krylov subspace methods (e.g., Conjugate Gradient, GMRES) [27,24]. Krylov methods approximate the solution within 17 1.5. Methods for Solving Linear Systems a subspace generated by the matrix and the initial guess, while multigrid methods leverage a hierarchy of discretizations to address multiple scales of behavior. While multigrid methods can be highly effective, they often require application-specific tuning and can face challenges when applied to highly anisotropic boundary layers. On the other hand, stationary methods work by finding the fixed point of a predefined mapping. 1.5.1 Direct Methods Direct methods solve a linear system in a finite number of steps. Despite their high computational cost and memory demands, they can sometimes outperform iterative solvers for small problems [29]. Direct methods are generally classified into elimination methods, such as Gaussian elimination and factorization methods, including LU and Cholesky factorizations. LU decomposition factorizes the linear system matrix A into the product of a lower triangular matrix L and an upper triangular matrix U , provided this decomposition exists. The system Ax = b can then be solved in two steps: Ly = b followed by Ux = y . The Cholesky decomposition is a particular case of LU decomposition for symmetric and positive-definite matrices, factorizing A as A = LLT , where L is a lower triangular matrix. Solving Ax = b with Cholesky requires solving Ly = b and then LTx = y . Cholesky decomposition is generally more efficient than LU decomposition, requiring fewer arithmetic operations and less memory [29]. 1.5.2 Multigrid Methods Multigrid methods operate on a hierarchy of grids, refining errors on different scales to accelerate convergence. On fine grids, iterative solvers like GaussSeidel or Jacobi are used to reduce high-frequency errors (smoothing), while low-frequency errors are addressed by transferring the problem to coarser grids. Grid transfer operations include restriction (coarse-grid representation) and prolongation (fine-grid interpolation). These methods follow structured cycles, such as the V-cycle or W-cycle, balancing computational efficiency and error reduction. However, their performance can degrade for highly anisotropic problems, often requiring specialized smoothers [3,31]. 1.5.3 Iterative Methods Iterative methods require less computational resources and memory than direct methods but can suffer from slow convergence or divergence when dealing with ill-conditioned matrices. Preconditioning strategies are frequently applied to mitigate these issues and improve solver efficiency [2,24]. Given a linear system Ax =b, its residual r(x)is defined as: r(x) := b−Ax. 18 1.5. Methods for Solving Linear Systems Most iterative methods attempt to approach the exact solution ¯ x by generating a sequence of approximations that gradually reduce r(x) to zero. However, due to floating-point errors, reaching the exact solution is often impractical, and a predefined tolerance level is usually employed to determine when to stop the iterations. The efficiency of iterative solvers depends on the properties of the system matrix and the choice of preconditioners, which will be explored in the following sections. Krylov Methods Krylov methods construct a sequence of approximations within a Krylov subspace: Kn(r0) = span{r0, Ar0, A2r0, . . . , An−1r0}, where r0 = Ax0−b is the initial residual. These methods, including Conjugate Gradient (CG) and GMRES, minimize the residual over this subspace, often achieving acceptable accuracy in far fewer iterations than the system size N . CG is one of the most effective solvers for symmetric and positive-definite matrices, though its efficiency strongly depends on preconditioning [2,24]. In computational fluid dynamics, highly anisotropic meshes induce stiffness, leading to ill-conditioned matrices. This stiffness slows down convergence significantly, necessitating robust preconditioners. Conjugate Gradient Method The CG method solves Ax = b for symmetric positive-definite A by minimizing the quadratic function: f(x) = 1 2xTAx −xTb, (1.37) which has a global minimum at the exact solution. Instead of simple steepest descent, CG uses mutually independent (A-conjugate) search directions, reducing the number of necessary iterations. The CG algorithm proceeds as follows: 19 1.5. Methods for Solving Linear Systems Algorithm 1 Conjugate Gradient Method 1: r=b−Ax0 2: d=r 3: δnew =rTr 4: δ0=δnew 5: while i<imax and ||r||2> ϵ2do 6: q=Ad 7: α=δnew/(dTq) 8: x=x+αd 9: if imod 50 = 0 then 10: r=b−Ax 11: else 12: r=r−αq 13: δold =δnew 14: δnew =rTr 15: β=δnew/δold 16: d=r+βd 17: i=i+ 1 This method efficiently converges in at most N iterations but often achieves a satisfactory solution in significantly fewer steps. The output criterion ||r||2> ϵ2 is used to determine when to stop the iterations, where ϵ is a predefined tolerance level. Another common choice is to stop when the relative residual is below a certain threshold, e.g., δnew/||b||< ϵ. Preconditioned Conjugate Gradient Method When the system matrix A is ill-conditioned, as often occurs in the discretization of the Laplacian operator on highly anisotropic meshes, the numerical stiffness of the problem increases, significantly degrading the convergence rate of iterative solvers. Preconditioning techniques can be applied to mitigate this issue. A preconditioner is a matrix that modifies the original linear system to improve numerical stability and convergence speed. It achieves this by approximating the system matrix in a way that reduces its condition number and improves eigenvalue clustering, thus making iterative solvers more efficient [26,16]. There are several types of preconditioners, classified based on their construction and application [2,24]: • Left preconditioning: This technique modifies the operator applied to x , solving the system M−1Ax = M−1b . This transformation improves spectral properties but requires careful treatment of convergence criteria. The trivial choice M = A yields a condition number of one and convergence in a single iteration, but solving M−1Ax would be as costly as solving Ax =bdirectly. 20 1.5. Methods for Solving Linear Systems • Right preconditioning: This technique transforms the solution vector directly by solving AM−1y = b and then recovering x via Mx = y . This approach ensures that the transformed system retains the same eigenvalues as the original. • Mixed preconditioning: This technique involves a two-sided transformation M−1 LAM−1 Ry = M−1 Lb , followed by solving MRx = y . It combines elements of both left and right preconditioning to achieve better conditioning. The choice of the most suitable preconditioning technique is highly problemdependent. For non-symmetric matrices, methods like GMRES are more appropriate [8,25,24], whereas, for symmetric matrices, the Conjugate Gradient (CG) method is often the best choice. Preconditioning must preserve symmetry for the CG method, which requires a symmetric system matrix. The naive approach of applying M−1A does not necessarily maintain symmetry, even if M and A are symmetric. Instead, we introduce a symmetric factorization. Since M is symmetric and positive-definite, there exists a matrix E such that M = EET . Using this, we transform the original system: E−1AE−Tˆx=E−1b, ˆx=ETx. The transformed system retains symmetry and positive definiteness, allowing it to be solved efficiently with the CG method. Notably, explicit computation of E is unnecessary; it suffices to efficiently compute the action of M−1 on a vector. The performance of a preconditioner M depends on how well it reduces the condition number of M−1A and clusters its eigenvalues. The challenge is selecting a preconditioner that balances effectiveness and computational efficiency. Thus, the Preconditioned Conjugate Gradient (PCG) method involves the following iterative steps: 21 1.5. Methods for Solving Linear Systems Algorithm 2 Preconditioned Conjugate Gradient (PCG) Algorithm 1: i= 0 2: r=b−Ax0 3: d=M−1r 4: δnew =rTr 5: δ0=δnew 6: while i<imax and ||r||2> ϵ2do 7: q=Ad 8: α=δnew/(dTq) 9: x=x+αd 10: if imod 50 = 0 then 11: r=b−Ax 12: else 13: r=r−αq 14: s=M−1r 15: δold =δnew 16: δnew =rTs 17: β=δnew/δold 18: d=s+βd 19: i=i+ 1 1.5.4 Different types of preconditioners The simplest preconditioner is the Jacobi preconditioner, also known as diagonal preconditioning. It uses only the diagonal entries of the system matrix A to construct a diagonal matrix M , which is then used as the preconditioner. This method, which effectively scales the quadratic form (1.37) along the coordinate axes, is computationally inexpensive and easy to implement, as the inversion of a diagonal matrix is trivial. However, its effectiveness is often limited to well-conditioned or diagonally dominant matrices, making it less suitable for systems with strong anisotropies or large condition numbers. A more advanced class of preconditioners relies on incomplete factorizations, such as Incomplete Cholesky (IC) and Incomplete LU (ILU). The Incomplete Cholesky preconditioner is used for symmetric positive-definite matrices and approximates A by factoring it into LLT , where L is a lower triangular matrix. Unlike full Cholesky factorization, IC limits fill-in to preserve the sparsity pattern of A , reducing storage and computational costs. The preconditioner is applied indirectly by solving LLTw = z using back-substitution. While IC is often effective, it can suffer from numerical instabilities, particularly for ill-conditioned or indefinite matrices. Similarly, LU decomposition and its incomplete variant, Incomplete LU (ILU), are widely used preconditioners, particularly for non-symmetric matrices. LU decomposition factors A into the product of a lower triangular matrix L and an upper triangular matrix U . However, direct LU factorization introduces significant fill-in, increasing memory and computational costs. The ILU method alleviates this by allowing only a controlled amount of fill-ins, producing an approximation A≈LU . While solving LUx = b is efficient, ILU does not 22 1.6. Linelet preconditioner directly yield the solution to Ax = b . Instead, the preconditioner M = LU is used within an iterative solver such as the conjugate gradient method. Despite their widespread use, IC and ILU may require tuning parameters, such as drop tolerance, to control sparsity and maintain numerical stability. Multigrid methods, another prominent class of preconditioners, address error components at different spatial resolutions by employing a hierarchy of coarser grids to efficiently eliminate low-frequency errors, while finer grids refine high-frequency components. These highly scalable methods can achieve near-optimal convergence rates for elliptic and parabolic problems, making them a widely used approach in CFD applications. However, their implementation can be particularly challenging for unstructured grids or highly anisotropic domains, where convergence may be significantly degraded or even lead to solver divergence. Algebraic Multigrid (AMG), a variant of multigrid that builds coarse grids algebraically rather than geometrically, is also widely used in CFD due to its adaptability to complex geometries and unstructured meshes. Other advanced preconditioners include Schwarz methods, which partition the domain into overlapping or non-overlapping subdomains and solve smaller subproblems locally to construct a global preconditioner. Sparse Approximate Inverse (SPAI) preconditioners are a class of preconditioners that approximate the inverse of the system matrix A by constructing a sparse matrix M , i.e., M≈A−1 , that retains the sparsity pattern of A . This allows for efficient computation while improving the convergence rate of iterative solvers like Conjugate Gradient (CG). However, SPAI preconditioners may not always provide significant improvements if the original matrix A is ill-conditioned. In such cases, additional techniques may be needed to improve the stability and efficiency of the preconditioner. It is also worth mentioning that the choice of the preconditioner is highly problem-dependent. A preconditioner that performs well in one case might perform poorly in another, whereas a different preconditioner could be more effective. Also, a suitable preconditioner must balance computational cost, memory requirements, and convergence performance to ensure computational efficiency and good scalability. Therefore, selecting an appropriate preconditioner for iterative solvers is as critical as choosing the solver itself. 1.6 Linelet preconditioner Among the wide variety of such preconditioners, the preferred option to work with in cases of high mesh anisotropy is the linelet preconditioner [13,18,19]. Another kind of linelet preconditioner is the streamline linelet preconditioner, developed in previous works [8]. Nevertheless, we find the former more suitable for the present work. In this method, the preconditioning matrix M is built by filtering out some of the coefficients of matrix A . Briefly, the main idea 23 1.6. Linelet preconditioner (a) Sketch for different linelets built in an unstructured mesh. (b) Tridiagonal structure of the linelet preconditioning matrix for the case shown in Figure 1.2a. Figure 1.2: Linelet structure derived from the strongest couplings in the mesh and the tridiagonal matrix obtained by exclusively considering these dominant connections. consists of building lines following the direction of the strongest couplings: the so-called linelets. These structures will only include nodes whose couplings are stronger than a certain tolerance. Also, a node can be part of at most one linelet, depending on the strength of its couplings: if the couplings are strong enough, the node will be included in a linelet; otherwise, it will not, but the same node will never be part of two or more linelets. Then, the preconditioning matrix is built by assembling the diagonal entries of the system matrix ( Aii ) and the non-diagonal entries ( Aij , i=j ) of the edges belonging to these lines of high couplings, any other coupling not included in any linelet is discarded. This means that for the mesh nodes that do not belong to any linelet, only the diagonal entry is preserved in the preconditioning matrix. In other words, the diagonal preconditioner is applied [14] for nodes not included in any linelet. Finally, if these points are renumbered following the linelets, the preconditioning matrix associated is block-tridiagonal. This method is particularly useful when ill-conditioning comes from a high anisotropy in the mesh. Schematically, the linelet preconditioning matrix is defined as follows: Mii =Aii i.e., Mhas the same diagonal as A Mij =Aij i=jif i, j belongs to a linelet Mij = 0 i=jif i, j does not belong to a linelet. (1.38) After a nodal renumeration following the linelets, the structure of M consists of sets of decoupled tridiagonal matrices, each associated with a linelet. As stated above, only the diagonal entry remains for nodes not associated with any linelet. Figure 1.2a shows an example grid where the strongest couplings are shown in red, exposing the linelet structure of couplings. The structure of the corresponding preconditioner Mis shown in Figure 1.2b. Given that the set of preconditioning equations can be written as s=M−1r, (1.39) 24 1.6. Linelet preconditioner where r is the residual and M the preconditioning matrix whose tridiagonal structure is shown in Figure 1.2b, we need an efficient way to compute M−1r . The set of decoupled tridiagonal equations in M allows a fast solution of the preconditioning equations in PCG by a direct method, where one tridiagonal system needs to be solved for each linelet. Hence, the Tridiagonal Matrix Algorithm (TDMA) (also known as the Thomas Algorithm) is suitable for this step. It is important to note that the linelet preconditioner can be constructed either algebraically by directly inspecting the linear system matrix A , as previously mentioned, or geometrically. In the algebraic approach, linelets are formed by following the direction of the strongest couplings in matrix A . Conversely, in the geometric approach, linelets are constructed by connecting the closest nodes. This is analogous to the procedure described in (1.38) , but instead of analyzing the linear system matrix A , it relies on inspecting the matrix defined as Dij =(1 ||xi−xj||2,i=j 0,i=j(1.40) where x i and x j represent the coordinates of nodes i and j , respectively. Here, Dij corresponds to the inverse of the distance between nodes i and j . Once the linelets are constructed, the entries used to form the preconditioning matrix M are extracted from the linear system matrix A. 1.6.1 Caracterization of the linelet preconditioner In our experience, the linelet preconditioner, designed to exploit the mesh anisotropy in its favor, has demonstrated improved performance with increasing aspect ratios. Specifically, a higher aspect ratio correlated with a reduced number of iterations required for convergence, indicating faster convergence rates (we define the aspect ratio as the geometrical aspect ratio, representing the ratio of the maximum to minimum length of a given element). Consequently, the presented results are particularly advantageous for scenarios characterized by higher aspect ratios. Before analyzing the convergence of the proposed method in detail, let us first examine the effect of aspect ratio variation on the convergence of the linelet preconditioner. We consider two meshes: first, a hybrid mesh with 0.49 million elements, where 49% are within the boundary layer, and second, a structured mesh with 0.5 million elements, where 100% are within the boundary layer. In this context, we varied the maximum aspect ratio of the boundary layer, maintaining a progression of 1.01, and compared the number of iterations the solver requires to reach a relative residual of 10 −6 . To isolate the effect of the maximum aspect ratio on convergence, we ran the case sequentially and without domain partitioning; in this case, our new approach to the linelet preconditioner, which we call Global Linelet Preconditioner (GLP), is equivalent to previous ones Local Linelet Preconditioner (LLP), where no communication scheme is implemented. Additionally, we compared its behavior with the diagonal preconditioner. 25 1.6. Linelet preconditioner (a) Case tun on a mesh with 0.48M elements, 49% of them within the boundary layer. The number of iterations for the linelet preconditioner ranges from 310 for an aspect ratio of 14 to 321 for an aspect ratio of 1024. (b) Case run on a mesh with 0.5M elements, all within the boundary layer. The number of iterations for the linelet preconditioner goes from 94 for an aspect ratio of 14 to 11 for an aspect ratio of 1024. Figure 1.3: Comparison of the number of iterations required to reach a tolerance of 10 −6 for the linelet and diagonal preconditioner in the PCG Method on meshes with varying maximum aspect ratio and boundary layer coverage. 26 1.7. Auxiliary methods 1.7 Auxiliary methods This section introduces the TDMA and the Schur-complement method, which will be used in the following sections to solve the preconditioning equations. In particular, the Schur-complement method will be used to parallelize the TDMA. 1.7.1 TDMA: Tridiagonal Matrix Algorithm This is a simple, straightforward method used to find a direct solution for any linear system whose matrix is tridiagonal. In this work, it will be used to solve equations (1.47) and (1.48). Suppose we need to solve the following problem:         b1c10 0 0 . . . 0 a2b2c20 0 . . . 0 0a3b3c30. . . 0 0 0 a4b4c4. . . 0 . . . . . . . . . . . . . . . . . . . . . 00000. . . bn                 x1 x2 x3 x4 . . . xn         =         b1 b2 b3 b4 . . . bn         . If there is no need to preserve the values of bi and di , then the following iteration solves analytically for xi: For i=2,3,...,n do: wi=ai bi−1 bi=bi−wici−1 di=di−widi−1. (1.41) Followed by the back substitution: xn=dn bn xi=di−cixi+1 bi , (1.42) where the last line in (1.42) is iterated over i=n−1, n −2, ..., 1. 1.7.2 The Schur-complement method The Schur-complement method simplifies solving linear systems by partitioning matrices into blocks and eliminating variables, making it especially effective for sparse matrices with block structures. The main idea behind this method is to decompose the linear system into smaller subsystems involving only a subset of variables. Once the reduced system is solved, the remaining variables are determined using back-substitution. This method is particularly useful for parallel computing, as it allows each processor to solve its subdomain plus an 33 1.7. Auxiliary methods interface problem with just one global communication. Suppose a linear system Ax = b , where x is distributed among P processors. Interface and inner nodes are defined so that the inner nodes of a given subdomain are only coupled to inner nodes of the same subdomain and to interface nodes. Interface nodes, on the contrary, couple inner nodes of different subdomains. Note that inner nodes from different subdomains are not coupled within this approach. We take it apart to form a new subset xs so that x is partitioned into P subsets plus an interface s. Thus, vector xcan be reordered as x= (x0, x1, ..., xP−1, xs)t.(1.43) It is important to note that by definition, given two subsets xp1 and xp2 , the system matrix does not link elements of these subsets together, but only between xpi and the interface terms xs , hence equation Ax = b can be rewritten as:       A0,00. . . A0,s 0A1,1. . . A1,s . . . . . . . . . . . . 0. . . AP−1,P −1AP−1,s As,0As,1. . . As,s             x0 x1 . . . xP−1 xs       =      b0 b1 . . . bP−1 bs       (1.44) By block Gaussian elimination, the solution to equation (1.44) is given by solving the following set of linear equations: ˜ bs=bs− P−1 X p=0 As,pA−1 p,pbp(1.45) ˜ As,s =As,s − P−1 X p=0 As,pA−1 p,pAp,s.(1.46) Finally, the solution to xis found as: ˜ As,sxs=˜ bs(1.47) and ˜ Ap,pxp=bp−Ap,sxs.(1.48) Therefore, once the linear system is reordered as (1.44) , interface nodes can be solved by applying equations (1.45) , (1.46) and (1.47) . Using the solutions for these interface nodes xs , equation (1.48) solves the system for the remaining nodes. Note that if matrix A does not change during the execution, then matrices As,s , Ap,p , and Ap,s can be computed only once at the preprocessing stage. Also, the product As,pA−1 p,p can be computed only once. The way to do this will 34 1.8. Conclusions be explained in detail in section 1.7.2. 1.8 Conclusions This chapter has introduced the motivations and challenges underlying the development of the linelet preconditioner, particularly in the context of high-fidelity CFD simulations. One of the primary limitations of conventional solvers lies in the treatment of highly anisotropic meshes, where traditional preconditioners struggle to maintain efficiency and convergence. The linelet preconditioner was developed to address this issue by leveraging the dominant coupling directions in the system matrix. However, despite its effectiveness in handling these problems, extending the method to parallel computing environments without degrading the solver’s convergence presents new challenges, particularly in handling inter-subdomain couplings. It is common practice to constrain linelets to a single subdomain, preventing them from spanning multiple partitions. While this simplifies the application of the preconditioner, it also results in the exclusion of inter-subdomain couplings from the preconditioning matrix, leading to a degradation in convergence. This effect becomes especially pronounced in highly anisotropic meshes with many domain partitions. As shown in this chapter, the loss of these couplings can increase the number of solver iterations required to achieve a given level of accuracy. To overcome these limitations, we propose a novel parallelization strategy for the linelet preconditioner, which allows linelets to extend across multiple subdomains without degrading the solver’s convergence. This approach is based on the Schur complement method, which enables efficient handling of inter-subdomain dependencies while maintaining the sparsity and structure of the preconditioning matrix. By adapting existing algorithms to incorporate these couplings, we ensure that the preconditioner remains effective regardless of the domain partitioning strategy used. This mitigates the degradation in convergence observed in traditional methods and enhances solver robustness and scalability in large-scale simulations. By allowing linelets to be distributed across multiple processes without sacrificing accuracy, we enable more efficient simulations of complex flow phenomena, particularly in cases where high Reynolds numbers and boundary layer effects play a critical role. The subsequent chapters will delve into the theoretical formulation, implementation details, and numerical experiments that validate the effectiveness of this approach. References [1] Richard M Beam and R.F Warming. ‘An implicit finite-difference algorithm for hyperbolic systems in conservation-law form’. In: Journal of Computational Physics 22.1 (1976), pp. 87–110. 35 1.8. Conclusions [2] Michele Benzi. ‘Preconditioning Techniques for Large Linear Systems: A Survey’. In: Journal of Computational Physics 182.2 (2002), pp. 418–477. issn: 0021-9991. doi:https://doi.org/10.1006/jcph.2002.7176.url: https://www.sciencedirect.com/science/article/pii/S0021999102971767. [3] William L. Briggs, Van Emden Henson and Steve F. McCormick. A Multigrid Tutorial (2nd Ed.) USA: Society for Industrial and Applied Mathematics, 2000. isbn: 0898714621. [4] W.R. Briley and H. McDonald. ‘Solution of the multidimensional compressible Navier-Stokes equations by a generalized implicit method’. In: Journal of Computational Physics 24.4 (1977), pp. 372–397. [5] J.A. Carlson, A. Jaffe, A. Wiles, Clay Mathematics Institute and American Mathematical Society. The Millennium Prize Problems. American Mathematical Society, 2006, pp. 57–67. isbn: 9780821836798. url:https: //books.google.es/books?id=7wJIPJ80RdUC. [6] Ramon Codina. ‘Pressure Stability in Fractional Step Finite Element Methods for Incompressible Flows’. In: Journal of Computational Physics 170.1 (2001), pp. 112–140. issn: 0021-9991. doi:https://doi.org/10.1006/ jcph.2001.6725.url:https://www.sciencedirect.com/science/article/pii/ S0021999101967257. [7] Ramon Codina. ‘Pressure Stability in Fractional Step Finite Element Methods for Incompressible Flows’. In: Journal of Computational Physics 170 (June 2001), pp. 112–140. doi:10.1006/jcph.2001.6725. [8] P. Córdoba Pañella. ‘Local preconditioning for parallel iterative solvers.’ PhD thesis. UPC, Facultat de Matemàtiques i Estadística, 2022. url: http://hdl.handle.net/2117/370230. [9] D. J. D. Mavriplis. ‘On convergence acceleration techniques for unstructured meshes’. In: 29th AIAA, Fluid Dynamics Conference.doi:10.2514/ 6.1998-2966. eprint: https://arc.aiaa.org/doi/pdf/10.2514/6.1998-2966. url:https://arc.aiaa.org/doi/abs/10.2514/6.1998-2966. [10] V.Raghvendra D.Bahuguna and B.V.R.Kumar. Topics in sobolev spaces and applications. Narosa Publishing House, 2002, p. 28. isbn: 1842650947. [11] Timothy A. Davis. Direct Methods for Sparse Linear Systems. USA: Society for Industrial and Applied Mathematics, 2006, p. 67. isbn: 0898716136. [12] O. Hassan, K. Morgan and J. Peraire. ‘An implicit finite-element method for high-speed flows’. In: International Journal for Numerical Methods in Engineering 32.1 (1991), pp. 183–205. [13] Rainald Löhner; Fernando Mut; Juan Raul Cebral; Romain Aubry; Guillaume Houzeaux. ‘Deflated preconditioned conjugate gradient solvers for the pressure-Poisson equation: Extensions and improvements’. In: International Journal for Numerical Methods in Engineering 2010-jun 22 vol. 87 iss. 1-5 87 (1-5 June 2010). doi:10.1002/nme.2932.url: libgen.li/file.php?md5=385da88d3a20687c5af11860062f9112. [14] Roux F. Houzeaux G. Mogoules F. Parallel scientific computing. Wiley Computer Engineering Series, 2015, pp. 59, 62, 65, 66, 86, 88, 98, 99, 100, 120. isbn: 9781118761687. 36 1.8. Conclusions [15] Thomas JR Hughes. The finite element method: linear static and dynamic finite element analysis. Courier Corporation, 2003, p. 118. [16] Noriyuki Kushida. ‘Condition Number Estimation of Preconditioned Matrices’. In: PloS one 10 (Mar. 2015), e0122331. doi:10.1371/journal. pone.0122331. [17] Ulrich Langer and Martin Neumüller. ‘Direct and Iterative Solvers’. In: Computational Acoustics. Ed. by Manfred Kaltenbacher. Cham: Springer International Publishing, 2018, pp. 66–67. isbn: 978-3-319-59038-7. doi: 10.1007/978-3-319-59038-7_5.url:https://doi.org/10.1007/978-3-31959038-7_5. [18] D. Martin and R. Löhner. ‘An implicit, linelet-based solver for incompressible flows’. In: AIAA 1992-668. 30th Aerospace Sciences Meeting and Exhibit. Jan. 1992. [19] DOROTHEE MARTIN and RAINALD LOEHNER. ‘An implicit, lineletbased solver for incompressible flows’. In: 30th Aerospace Sciences Meeting and Exhibit.doi:10.2514/6.1992-668. eprint: https://arc.aiaa.org/doi/pdf/ 10.2514/6.1992-668.url:https://arc.aiaa.org/doi/abs/10.2514/6.1992668. [20] D. J. Mavriplis. ‘Directional Agglomeration Multigrid Techniques for High-Reynolds-Number Viscous Flows’. In: AIAA Journal 37.10 (1999), pp. 1222–1230. [21] D. J. Mavriplis. ‘On the accuracy of the pseudocompressibility method in solving the incompressible Navier-Stokes equations’. In: Applied Mathematical Modelling 11.1 (1987), pp. 35–44. [22] J. Perot. ‘An Analysis of the Fractional Step Method’. In: Journal of Computational Physics 108 (Sept. 1993), pp. 51–58. doi:10.1006/jcph. 1993.1162. [23] Ravi Ramamurti and Rainald Löhner. ‘A parallel implicit incompressible flow solver using unstructured meshes’. In: Computers and Fluids 25 (Feb. 1996). doi:10.1016/0045-7930(95)00032-1. [24] Y. Saad. Iterative Methods for Sparse Linear Systems. 2nd. USA: Society for Industrial and Applied Mathematics, 2003, pp. 4, 5, 57, 63, 69, 75, 82, 83, 90, 92, 100, 149. isbn: 0898715342. [25] Youcef Saad and Martin H. Schultz. ‘GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems’. In: SIAM Journal on Scientific and Statistical Computing 7.3 (1986), pp. 856–869. doi:10.1137/0907058. eprint: https://doi.org/10.1137/0907058.url: https://doi.org/10.1137/0907058. [26] Yousef Saad. Iterative Methods for Sparse Linear Systems. Second. Society for Industrial and Applied Mathematics, 2003. doi:10.1137/ 1.9780898718003. eprint: https://epubs.siam.org/doi/pdf/10.1137/ 1.9780898718003.url:https://epubs.siam.org/doi/abs/10.1137/1. 9780898718003. [27] Jonathan R Shewchuk. An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. Tech. rep. USA, 1994. 37 1.8. Conclusions [28] Orlando Soto, Rainald Lohner and Fernando Camelli. ‘A Linelet preconditioner for incompressible flow solvers’. In: International Journal of Numerical Methods for Heat & Fluid Flow - INT J NUMER METHOD HEAT FL F 13 (Feb. 2003), pp. 133–147. doi:10.1108/ 09615530310456796. [29] JM Vargas-Felix and S Botello-Rionda. ‘Parallel direct solvers for finite element problems’. In: Comunicaciones del CIMAT, I-10-08 (CC) (2010). [30] Jacob Waltz. ‘Unstructured multigrid for time-dependent incompressible fluid flow. PhD Thesis, School of Computational Sciences of George Mason University’. In: (Jan. 2000). [31] P. Wesseling. Introduction to Multogrid Methods. Tech. rep. 1995. [32] P. Wesseling and C.W. Oosterlee. ‘Geometric multigrid with applications to computational fluid dynamics’. In: Journal of Computational and Applied Mathematics 128.1 (2001). Numerical Analysis 2000. Vol. VII: Partial Differential Equations, pp. 311–334. 38 CHAPTER 2 Efficient Parallelization of the Tridiagonal Matrix Algorithm Part of the contents of this chapter have been published as: R. de Olazábal, R. Borrell and O. Lehmkuhl. ‘An algebraic global linelet preconditioner for incompressible flow solvers’. In: Journal of Computational Physics 514 (2024), This method is primarily based on the previous work done by O. Soto, R. Löhner, and F. Camelli [2] on a linelet preconditioner for anisotropic meshes. In a way, it serves as an extension of it and its natural continuation. One of the main challenges of the linelet preconditioner is that it requires inverting a tridiagonal matrix whose entries may be distributed across many processes. Therefore, the usual approach is to discard interdomain couplings. In this work, we propose a method that allows for these couplings to be taken into account, thus improving the preconditioner’s performance. Two primary challenges must be addressed to achieve this goal. The first to perform a Tridiagonal Matrix Algorithm (TDMA) on the preconditioning matrix M when its entries are distributed across multiple processes. The second is efficiently assembling M while preserving the necessary couplings between subdomains. This chapter focuses on the first challenge: parallelizing the TDMA by employing the Schur Complement Method. This approach enables the efficient solution of the linear system (1.39) in parallel while maintaining the tridiagonal structure of the preconditioning matrix M . We also present a technique to efficiently compute each subdomain’s contribution to the r.h.s., minimizing redundant computations and communication overhead. By leveraging the inherent sparsity of the tridiagonal system and optimizing interprocess communication, we demonstrate how to significantly reduce computational costs while maintaining the accuracy and robustness of the preconditioner. Finally, note that the method presented in this chapter is not restricted to the linelet preconditioner. It can be applied to any tridiagonal matrix whose entries are distributed across multiple processes. However, throughout this 39 2.1. Various approaches chapter, we focus on linelets and the linelet preconditioner without loss of generality, as the connectivity graph of any tridiagonal matrix forms a 1D structure, which can be interpreted as a linelet. 2.1 Various approaches Linelet preconditioners have been widely used in the literature. Since parallelizing them when the entries of the preconditioning matrix are spread over many processes is not trivial, one of two approaches is taken when dealing with these types of preconditioners in parallel solvers. 2.1.1 Previous approaches The first approach is to let each subdomain define its linelets. In this approach, no process shares information with another regarding the linelets. In this case, the domain partition is independent of the linelets. The absence of communications implies that each subdomain may contain sections of different linelets, but these sections are treated as separated linelets on their own. This approach, as shown in figures 1.8-1.10, can limit the numerical performance of the preconditioner. The second one is to adapt the mesh partition to the size of the linelets within each subdomain. Since linelets can include more elements than in the previous case, and thus more couplings are taken into account, a major correspondence between matrices A and M is achieved. Nevertheless, adapting the mesh may lead to a poor load balance or, due to the geometry of the problem, it may not be possible to accommodate every linelet within a single subdomain. In this situation, it will still be necessary to cut some of the linelets, leaving us in a similar situation to the previous case. 2.1.2 Our approach: agnostic to the domain partition In this work, we try a different approach: the domain partition is defined independently of the linelet geometry, but a communication step is implemented within the preconditioning step. Therefore, the linelet structures, and therefore the preconditioner, are agnostic to the domain partition. This is done so that different sections of each linelet can share information. Hence, the number of iterations will be reduced compared to the case where these communications are not allowed (since, algebraically, it will be equivalent to solving the linear system shown in equation (1.39) sequentially, letting the linelets be as long as possible and achieving a major correspondence between matrices A and M ). Nevertheless, this does not necessarily imply a reduction in the total solver time since a communication step needs to be included, and thus, it will generate an overhead. Therefore, a trade-off needs to be considered between the reduction in the number of iterations and the increase in cost per iteration. For clarification purposes, we will redefine the terms “linelet” and “slice of a linelet” to be used in this work. As shown in Figure 2.1, we will use the term “slice” to refer to the local structure given by the nodes located within 40 2.2. Parallelization of the TDMA using the Schur-Complement Method Figure 2.1: Definition of the terms "slice" and "linelet": We introduce the term "slice" to refer to the local structure of nodes that are connected within a single subdomain. In contrast, "linelet" refers to the global structure arising from many slices’ connection. Interface nodes are depicted as empty nodes, and inner nodes as filled ones. one subdomain. “Linelet”, on the other hand, will refer to the global structure given by coupling different slices. Note that a slice can be composed of a single node, and a linelet can be composed of a single slice. If a whole linelet is located within a single subdomain, equation (1.39) could be directly solved via the TDMA previously explained. Nevertheless, this work aims to develop a method able to solve equation (1.39) even if a linelet is spread over many processes, i.e., even if entries of the tridiagonal matrix M shown in figure 1.2b are distributed over many processes. Therefore, we want a way to perform a TDMA even if the entries of the matrix are spread over many processes. The Schur Complement Method was found to be a suitable solution for dealing with this case, parallelizing the TDMA. 2.2 Parallelization of the TDMA using the Schur-Complement Method This section addresses the parallelization of the TDMA using the SchurComplement method. We start this discussion by focusing on efficiently computing each slice’s contribution to the r.h.s.. In particular, we avoid unnecessary computations and communication costs. By profiting from the tridiagonal structure of the matrix and careful handling of communication between subdomains, we demonstrate how to calculate the contributions without requiring a full TDMA operation at every iteration. Once the contributions are computed, the next step is to gather them in a central process and solve the system. The tridiagonal structure of the matrix allows us to apply the Schur-Complement method to solve the system, breaking the problem into smaller slices and distributing the workload across different processes. 2.2.1 An efficient way to compute each slice’s contribution to the r.h.s. Recall from Figures 1.2a and 1.2b that each block in the preconditioning matrix M corresponds to a linelet, with couplings between nodes in the linelets limited to first neighbors. By dividing a linelet into sections and correctly identifying 41 2.2. Parallelization of the TDMA using the Schur-Complement Method the interface nodes, we can apply the Schur Complement Method to solve the linear system (1.39) . The slicing of linelets is determined by the domain partition, with the interface nodes chosen as the first ones after the interdomain interface, as shown in Figure 2.1. In this figure, interface nodes are represented as empty circles, while inner nodes are depicted as filled circles. Importantly, inner nodes from different subdomains are not coupled. This section explains the algorithm used to compute the contributions to the linear system (1.39) when Mis tridiagonal, and its entries are distributed across multiple subdomains. From equations (1.45) and (1.46) , it follows that each process managing a slice of a linelet must compute its contribution to equation (1.47) . Specifically, each process needs to evaluate As,pA−1 p,pbp and As,pA−1 p,pAp,s . However, if the system matrix remains constant during the system’s evolution, As,pA−1 p,pAp,s can be computed once during preprocessing. In contrast, As,pA−1 p,pbp must be updated at every iteration. A direct approach would involve using a TDMA to compute A−1 p,pbp and then multiplying it by As,p , but this would require performing a TDMA computation every time ˜ bs is needed. Instead, this section demonstrates that no TDMA is necessary to compute ˆ bs from equation (1.45) . At most, only two columns of As,pA−1 p,p need to be computed at the start, after which two dot products with bp suffice to calculate As,pA−1 p,pbp during each iteration. As highlighted in the Schur algorithm (equations (1.45) – (1.48) ), when the matrix remains unchanged between iterations, the product As,pA−1 p,p , along with As,s , can be computed once at the preprocessing stage. Taking into account the structure of the Schur matrices Ap,p , As,p , and the tridiagonal form of the preconditioner matrix M (Figure 1.2b), this product results in a sparse matrix. At least one 2 × 2diagonal block is nonzero. It can be explicitly computed, such as by applying TDMA to extract the required components of As,pA−1 p,p , as referenced in equations (1.45) and (1.46). Direct computation of the product As,pA−1 p,p Consider the example of a linelet spread over three processes as shown in figure 2.2. Note that we may need a proper renumeration of the nodes to achieve a tridiagonal structure. Suppose nodes 1to 5belong to process 0, nodes 6 to 9belong to process 1, and nodes 10 to 12 belong to process 2. We divide this matrix into different blocks following the linelet partition imposed by the domain partition, and we reorder them as shown in equation (1.44) . The result can be seen in figure 2.2. In this case, we have four diagonal blocks, A0,0 , A1,1 , A2,2 , and the interfaces block As,s , together with six non diagonal blocks, As,0 , As,1, and As,2(and their respective transposes). Since matrix As,p has at most two nonzero entries, the vector resulting from operating As,pA−1 p,pbp , will also have at most two nonzero elements. Hence, instead of computing all of them, we can take advantage of the fact that we already know the structure of As,p . Thus, we can compute only the nonzero entries while saving computational resources. A straightforward approach to compute As,pA−1 p,pbp would be to apply a TDMA to obtain A−1 p,pbp , and then multiply it with As,p . However, this would 42 2.3. Conclusions 2.3 Conclusions This chapter addresses the proposed method to parallelize the TDMA when applied to matrices whose entries are distributed across multiple processes. We developed a strategy to preserve interdomain couplings while maintaining computational efficiency. This approach offers a solution agnostic to domain partitioning, improving the correspondence between the system matrix A and the preconditioner M , reducing the number of iterations required for convergence. The work introduces an algebraic and parallel framework for implementing a Linelet Preconditioner to solve Poisson’s equation, generalizing the conventional approach by allowing linelets to span multiple subdomains. This ensures a more balanced workload distribution across computational resources. The algebraic approach of this method allows for broader application to matrices with similar properties. The flexibility of the algebraic method also opens up opportunities for extending its use beyond the Poisson equation to any matrix exhibiting linelet structure, regardless of geometry. Additionally, the preconditioning matrix can also be assembled based on geometric considerations, which emphasizes the capability of this approach to be adapted to different problems, with the added benefit that the domain decomposition can be independently optimized for load balancing, thereby minimizing computational overhead. In order to increase efficiency, we devised a strategy to compute each subdomain’s contribution to the r.h.s. of the system while taking advantage of the sparsity of the tridiagonal matrix to minimize redundant computations and communication overhead. We then introduced a parallel scheme for solving the resulting system, mapping linelets to processes in a way that balances computational load while minimizing communication costs. It is worth mentioning that this technique is not limited to linelet-based preconditioners, as it can be applied to any tridiagonal system distributed across multiple processes. In the next chapter, we will extend this work by addressing the efficient assembly of the preconditioning matrix M , ensuring that the necessary couplings between subdomains are preserved while maintaining scalability in large-scale parallel simulations. References [1] Yousef Saad. Iterative Methods for Sparse Linear Systems. Second. Society for Industrial and Applied Mathematics, 2003. doi:10.1137/ 1.9780898718003. eprint: https://epubs.siam.org/doi/pdf/10.1137/ 1.9780898718003.url:https://epubs.siam.org/doi/abs/10.1137/1. 9780898718003. 49 2.3. Conclusions [2] Orlando Soto, Rainald Lohner and Fernando Camelli. ‘A Linelet preconditioner for incompressible flow solvers’. In: International Journal of Numerical Methods for Heat & Fluid Flow - INT J NUMER METHOD HEAT FL F 13 (Feb. 2003), pp. 133–147. doi:10.1108/ 09615530310456796. 50 CHAPTER 3 Assembling the preconditioning matrix Part of the contents of this chapter have been published as: R. de Olazábal, R. Borrell and O. Lehmkuhl. ‘An algebraic global linelet preconditioner for incompressible flow solvers’. In: Journal of Computational Physics 514 (2024), and disseminated in the following conferences: the 32nd Parallel Computational Fluid Dynamics Conference, Nice, France, 2021, the 15th JLESC Workshop, Bordeaux, France, 2023, the 1st Math 2 Product (M2P) conference for Emerging Technologies in Computational Science for Industry, Sustainability and Innovation, Taormina, Sicily, 2023, and the 23rd IACM Computational Fluids Conference, Santiago de Chile, 2025. In the previous chapter, we introduced an efficient parallelization strategy for the TDMA when applied to distributed matrices. This approach can be implemented in the Linelet Preconditioner, extending the conventional approach to allow linelets to span multiple subdomains. This generalization ensures a more balanced workload distribution and expands the method’s applicability to matrices with similar properties. Previous approaches to linelet preconditioners involve allowing each subdomain to define its own linelets independently. This can hinder numerical performance due to the lack of communication between subdomains or needing to adapt the mesh partition to the linelet geometry. While the latter improves matrix correspondence, it often results in poor load balancing and further challenges with subdomain compatibility. This chapter presents a method for constructing and combining the slices into larger linelet structures. Recall that we define a "slice" as the local structure within a subdomain, while a "linelet" refers to the global structure formed by 51 3.1. Assembling the slices merging multiple slices across subdomains. Special attention must be paid to connecting linelets from different domains, constructing the linear system for the linelet preconditioner when matrix entries are distributed across multiple subdomains, and efficiently performing the necessary inter-domain communications. 3.1 Assembling the slices The assembly of the preconditioning matrix M is composed of two steps. In the first step, each subdomain independently defines its set of slices. This is crucial to obtain a reasonable approximation of the inverse of the system matrix A , as the effectiveness of a linelet preconditioner depends strongly on the proper choice of the set of linelets. Several known methods can be used to perform this step. For example, in [1], this is done by following streamlines, while in [3], it is done by following the strongest couplings. We will stick with the latter in this work, although the same implementation can be used for the former. In the second step, linelet slices that belong to different subdomains are joined together. However, no method has been found in the literature that performs this step. Therefore, one of the primary contributions of our work is developing a method that can efficiently perform this operation, i.e., to extend the linelets beyond the boundary between different subdomains. 3.1.1 Algebraic slice assembly within each subdomain Let us start this section by briefly mentioning the method used to create the local slices. Briefly, this step starts by randomly selecting a node as the source point for a slice; then, it continues growing this slice until it can no longer include additional nodes. Then, another node, not yet part of any slice, is chosen as the source for a new slice. This process is repeated until growing additional slices is no longer possible. In the first step, every subdomain builds its own set of linelet slices independently. As mentioned above, this is done by keeping only the strongest couplings while discarding the rest (see conditions (1.38) ). The algorithm picks one node randomly and then grows the slice following the strongest couplings. This process is repeated until no more slices can be grown within that subdomain. Schematically, it consists of the following steps: 1. A node that does not already belong to any slice is picked at random. Let us call it node i. 2. Then, we look for the strongest coupling, i.e., the coupling Aij that satisfies |Aij| |Aim|≥αfor Nc−2values of m,m=j,i=j,α≥1.(3.1) If such node j does not belong to any slice, it is taken to be the second node of the slice. Note that the longest linelets are obtained by setting 52 3.1. Assembling the slices α = 1. This value can be adjusted to couple only those nodes with a sufficiently strong coupling. The higher the value of α , the stronger the couplings that the linelets admit. Within this step, it also should be checked that for node j , its coupling with i is among its two strongest couplings. The reason why the next node in the slice should satisfy this condition for Nc− 2values of m , m=j , i=j and not just for every value of m=jwill be explained in the following section. 3. Then we move to node j and repeat the process. It is important in this step to ignore the coupling with node i , which already belongs to the linelet. 4. This process is repeated until a) we come across a node already included in a linelet, b) condition (3.1) is not satisfied, or c) the next node is part of the halo (in Figure (3.1d) this process stops when node j is coupled to node kin the halo). 5. If one of the previous things happens, we go back to the first node i and start searching again but now ignoring the nodes already belonging to the linelet (Figure 3.1d, the search is started again from node i ignoring the previous coupling with node j). 6. The process is iterated until previous conditions are no longer satisfied. This process is also exemplified in Figure 3.1. When dealing with ill-conditioned matrices arising from highly stretched boundary layers, it is worth noticing that some of the previous steps could be avoided by setting to grow the linelets in the orthogonal direction to the mesh stretching. Nevertheless, we want this method to be algebraic and avoid any information regarding the geometry of the mesh. At the end of this step, each process will have all the information related to the structure and couplings of its own slices, together with the couplings between them and the halo. This last piece of data is important because it will allow us to expand the slices (i.e., the local section of the linelets) beyond the subdomain interface. Also, at this point, each subdomain has information regarding its own set of slices. However, they have no information on what happens in other subdomains nor how the communications must be performed. We need to define the communication scheme to let different slices share information. Remarks on choosing the next node in the slice Previously, we noted that when growing the slices, it is important to check that for the next no node j to be included in the slice, its coupling with the previous node i in the slice is among the two strongest couplings for both of them. Now, we give a more detailed explanation of why this condition is necessary and how to implement it. Let us start with some considerations about equation (3.1) . Suppose the case shown in Figure 3.2, where the first node picked at random is 53 3.1. Assembling the slices (a) A node not part of any existing linelet slice is selected randomly. Then, the strongest coupling needs to be identified. (b) We verify that for node j , its connection with i is among its two strongest couplings. (c) We move to node j and repeat the process. (d) As node j is linked to the halo, the process halts, restarting from node i but excluding the couplings already included in the slice. (e) The process is repeated until equation (3.1) is no longer satisfied. (f) The same process is repeated for every node in the subdomain. Figure 3.1: Iterative process to build slices (i.e., the local section of the linelets within each subdomain). This iterative process is performed for every node in the subdomain until conditions (3.1) are no longer satisfied, resulting in the construction of a complete set of slices. 54 3.1. Assembling the slices Figure 3.2: Sketch showing nodes and couplings for an example mesh. Nodes h , iand jshould belong to the same slice. node i . Then, to grow the slice, we need to find the strongest coupling. Suppose in this case the strongest coupling is Aij , between nodes i and j (i.e., Aij > Aim ∀m=j,i=j, with j). In the second step, we need to check whether coupling Aij satisfies equation (3.1). Suppose in this case we take α= 2, so condition (3.1) reads |Aij| |Aim|≥2,m=j. (3.2) Figure 3.2 shows that we should include nodes h , i , and j in the slice under construction. However, equation (3.2) will not (and should not, as we shall see) hold for the couplings between these three nodes. Since each of them has two strongest couplings, the previous condition will not be valid (for example, for node i , we have Aij/Aih = 10 / 9 < 2). Nevertheless, to have a slice joining these nodes, condition (3.2) must be changed to apply only for couplings not included in the slice. Since each node can be coupled to at most two other nodes, condition (3.2) should not hold for these two couplings; hence, (3.2) needs to hold for Nc− 2couplings, where Nc is the total number of nodes coupled to node i. Condition (3.1) then becomes |Aij| |Aim|≥αfor Nc−2values of m,m=j,i=j . (3.3) The previous statement is valid whether node h should be included in the current slice. After finding that node i has its strongest coupling with j , one might consider taking node j as the next node of the slice. Nevertheless, we need to be careful to avoid the case where the couplings between Ajk and Ajl were much stronger than the coupling Aij . In other words, it could be the case that node i has its strongest coupling with node j , but coupling with i is not among the two strongest couplings of node j ( Aij > Aim ∀m=j , i=j , but Aji ≤Ajm ∀m=i , j=i ). Hence, we could have the cases shown in Figures 55 3.1. Assembling the slices (a) Sketch of a case where, by building the slices following the strongest couplings, node j should belong to the same linelet as nodes hand i. (b) In this case, node j should not be included in the same linelet as nodes h and i . However, from the perspective of node i , this situation is indistinguishable from the one depicted in Figure 3.3a. Figure 3.3: Illustration highlighting the impact of strongest couplings on linelet construction. In one scenario, node j is expected to belong to the same linelet as nodes h and i (Figure 3.3a). However, in a different scenario, node j should not be included in the same linelet as nodes h and i (Figure 3.3b), even if node i has the same couplings as before. In other words, whether or not node j should be part of this slice is determined by the couplings of node j rather than those of node i. 3.3a and 3.3b. For example, from these figures, it can be seen that even if the coupling between nodes i and j is dominant when taking into account all nodes coupled to i , for node j to be able to be included in the same slice, the same coupling between i and j needs also to be dominant when restricted to all nodes coupled to j . Although in both figures, node i “sees” the same couplings, in Figure 3.3a nodes h , i and j should be included in the same slice, but in Figure 3.3b one slice should contain nodes h and i , and another one should contain nodes k , j and l . The only difference in both scenarios is how node j (and not i ) is coupled to its neighbors. This is why, before adding node j node to the slice, we must check that from the point of view of j , coupling Aij is also among the two dominant couplings. In other words coupling Aij will be included in the slice if: 1. Aij is among the two strongest couplings of node i : |Aij | |Aim|≥α for Nc− 2 values of m,m=j,i=j. 2. Aij is among the two greatest couplings of node j:|Aij | |Ajm|≥αfor Nc−2 56 3.1. Assembling the slices values of m,m=i,i=j. Algorithm 3 outlines the algorithm used to build the slices. Here, we denote αtol the tolerance imposed in condition (3.3) . The main idea behind this approach to slice assembly is to explore the nodes and their connectivity given by matrix A and to identify and assemble slices based on the relative weight of the couplings Aij. Algorithm 3 Build slices 1: for i in nodes do 2: if node i is not included in any slice then 3: IsAvailable=False 4: Look for the strongest coupling Aij .▷Let jbe the node with the strongest coupling. 5: Compute the minimum value of αfor |Aij | |Aim|≥αfor Nc−2values of m,m=j,i=j. 6: if (Aij is not in the halo) and (is among the two strongest coupling of node j )then 7: IsAvailable=True 8: SecondNode ←j 9: else 10: Look for the second strongest coupling Aij .▷Let jbe this node. 11: Compute the minimum value of αfor |Aij | |Aim|≥αfor Nc−2values of m,m=j,i=j. 12: if Aji is among the two strongest coupling of node jthen 13: IsAvailable=True 14: while (IsAvailable) and (α≥αtol) and (node jis not included in any slice) do 15: Add node jto the slice 16: Look for the strongest coupling Ajk.▷Let kbe this node. 17: Compute the minimum value of αfor |Ajk | |Ajm|≥αfor Nc−2values of m,m=j,j=k . 18: if (node kis in halo) and (α≥αtol)then 19: The coupling between node jand the halo is stored. 20: if Ajk is among the two strongest coupling of node kthen 21: IsAvailable=True 22: j←k 23: if (kis in the halo) or (α < αtol) or (IsAvailable==False) then 24: j←SecondNode This algorithm begins by iterating through each node (line 1). Letting i denote the starting node, the process checks whether i is not part of an existing slice (line 2). If not, i is designated as a potential starting point. Then, in line 3, an auxiliary boolean IsAvailable is set to False. This variable will determine whether the growing process can be repeated from the next node in the slice. 57 3.1. Assembling the slices This variable will be set to True if the next node j is not in the halo and the coupling Aji is among its two strongest couplings. This is because if j is in the halo, the algorithm cannot continue growing the slice. Also, if j is not in the halo, but the coupling Aji is not among its two strongest couplings, the algorithm cannot continue growing the slice either. Subsequently, the algorithm identifies the strongest coupling from i , denoted as node j (line 4). If j is not in the halo, the algorithm verifies that the coupling Aij satisfies condition I) (3.3) (line 5). If satisfied, the algorithm verifies whether the coupling Aji is among the two strongest couplings of node j (lines 6-8, remember that since we assume symmetry, then Aij = Aji ). On the contrary, if node j is in the halo or the coupling conditions are not satisfied, the algorithm restarts the search from node i (lines 10-14). Then, the algorithm enters a loop to extend the slice, adding nodes based on the strongest coupling from the current node j (lines 16-30). The length of this slice and the nodes included in it are updated at every iteration (line 25). The loop continues until a) the current node is already part of a slice, b) the coupling conditions are no longer satisfied, or c) the next node is in the halo. If a termination condition is met, the algorithm returns to the initial node i and restarts the search in the other direction (line 28), ignoring nodes already included in the slice. The process iterates across all nodes (line 1), repeating the process until no node i is left unexplored. The process terminates when all nodes have been addressed, and every slice is built. At the end of this step, each process will have all the information related to the structure and couplings of its slices, together with the couplings between them and the halo. This last piece of data is crucial because it will allow us to expand the slices (i.e., the local section of the linelets) beyond the subdomain interface. These couplings are indicated in Figure 3.1f as dotted red lines. Also, at this point, each subdomain contains information about its slices. Still, it lacks information on the slices other subdomains have built or the required communication procedures (i.e., how the communications must be performed). The next step is to define the communication scheme. 3.1.2 Geometric slice assembly within each subdomain In the previous section, we discussed how to build the slices within each subdomain by inspecting the system matrix A . Other methods can be used to build the slices, such as following streamlines [1], using the strongest couplings, or using a geometric approach. Suppose we follow a geometric approach instead of preferring an algebraic one. In that case, we can assemble the slices by computing the distance between nodes and selecting the closet node (note that it is still important that the equation (3.1) holds for any new node to the slice) or, on the other hand, we can define a matrix of distances between nodes and use this matrix to build the slices. In our case, we generate linelets via a geometric approach by defining a matrix D as in (1.40) . Then, we follow the same steps as in the previous section, but matrix D is now used to build the slice instead of the linear system matrix. In this case, matrix D is only used to build the slices, but the preconditioning matrix M entries are still filled with the corresponding matrix A . In other words, the corresponding entries of the Schur-Complement matrices ( (1.45) - (1.48) ) are filled with the entries of the 58 3.2. Halo cleaning Algorithm 5 Algorithm used to symmetrize interdomain couplings. 1: Buffer ←IdGlobal 2: Halo update on Buffer 3: BufferHaloUpdated ←Buffer 4: for i in haloCouplings do 5: Couplings(i) ▷Looks for the couplings to the halo 6: Buffer[i] ←buffer[FirstNodeCoupling] 7: Halo update on Buffer 8: for i in haloCouplings do 9: Couplings(i) 10: if IdGlobal[i]==Buffer[FirstNodeCoupling] then 11: This coupling remains 12: else if nCouplings(i)==2 then 13: if IdGlobal[i]==Buffer[LastNodeCoupling] then 14: This coupling remains 15: for i in haloCouplings do 16: Couplings(i) 17: BufferHaloUpdated[i] ←BufferHaloUpdated[LastNodeCoupling] 18: Halo update on BufferHaloUpdated 19: for i in haloCouplings do 20: Couplings(i) 21: if IdGlobal[i]==BufferHaloUpdated[FirstNodeCoupling] then 22: This coupling remains 23: else if nCouplings(i)==2 then 24: if IdGlobal[i]==BufferHaloUpdated[LastNodeCoupling] then 25: This coupling remains 26: Every other coupling is discarded through all nodes, which at the same time are part of a slice and are coupled to the halo. For each one of these nodes, line 5 searches for the node in the halo to which it is coupled. Line 6, the entry associated with that node in Buffer is replaced by the value of the node to which it is coupled, which we call FirstNodeCoupling. In Figure 3.7, this is performed in the step denoted by Buffer[i] ← buffer[FirstNodeCoupling]. This is what happens, for example, for node i in the left column where, since it is coupled to node m (see Figure 3.7), its ID is replaced by the received ID of node m . Suppose it is coupled to two other nodes in the halo - such as the case with node l in Figure 3.5 - only one of these couplings is considered in this step. This is depicted in the left column of Figure 3.7, where the ID of node l -which is coupled to nodes q and r - is replaced by q (the idea of having a second array called BufferHaloUpdated is to repeat this process with those couplings that have not been considered previously). Then, a second halo update is performed, and the process of replacing entries in Buffer with the values in the halo to whom they are coupled is repeated. For all those nodes whose own ID has been recovered, FirstNodeCoupling should not be discarded. This is shown in Figure 3.7 after the second Buffer[i] ← buffer[FirstNodeCoupling]: when the received ID (and color) matches the one in array IdGlobal, then FirstNodeCoupling 65 3.3. Growing the linelets: joining slices from different subdomains Figure 3.8: Final stage of the case shown in Figure 3.5. Only the symmetric interdomain couplings remain, while all the non-symmetric couplings were discarded after applying the halo cleaning process described in this section. is not discarded. Then, the same steps performed between lines 4 and 18 are again performed between lines 19 and 33, with the only change that if a node is coupled twice to the halo, then LastNodeCoupling is used instead of FirstNodeCoupling. Also, this is performed with the array BufferHaloUpdated defined in line 21 (after the first halo update in Figure 5). Once more, for every node whose ID has been recovered, LastNodeCoupling has not been discarded. On the contrary, FirstNodeCoupling and LastNodeCoupling are discarded for every coupling whose ID has not been recovered. For example, in Figure 5, since nodes i , o do not recover their ID at final stages 1 and 2, their coupling to the halo must be discarded. This transforms the initial state shown in Figure 3.5 to the one shown in 3.8. At this point, every coupling between nodes from different subdomains is symmetric. This will be crucial for the next step, which consists of growing the linelets by joining slices from different subdomains. This is explained in the next section. 3.3 Growing the linelets: joining slices from different subdomains In this section, we will describe how to couple slices from different subdomains into larger structers called linelets in the context of this work. This is one of the most important steps in the construction of the preconditioner, as it is the step where the communication between different subdomains is defined. Also, the operations performed here will be the reason why this preconditioner is agnostic to the domain partition. There are two main ways to couple slices from different subdomains. The first one is to select a specific set of slices to be the starting point of their respective linelets. One reason to do this is if we want the linelets to start at the boundary of the mesh. In such a case, we would select the slices containing the nodes at the boundary of the mesh to be the starting point of their respective linelets. The second way is to let the slices join one another without specifying any starting point. This is the approach we will follow in 66 3.3. Growing the linelets: joining slices from different subdomains this work. Nevertheless, for the sake of completeness, we will explain both methods. 3.3.1 Selecting a specific set of slices to be the starting point of their respective linelets In this approach, the linelets are grown from a predetermined set of slices. We can take one of the inputs of the solver to be the ID of the nodes candidates to belong to the first slice of their linelet (considering that they do not need to be at one extreme of their slice, they may be, for example, at the middle of a slice), any slice containing such nodes will be considered the first one of their linelet. The next step consists of the coupling of slices from different subdomains. This is done in an iterative process by “growing” the linelets starting from certain slices. Since one of the inputs was the ID of the nodes candidates to belong to the first slice of their linelet (considering that they do not need to be at the end of their slice), any slice containing such nodes will be considered the first one of their linelet. This step can be better understood by considering the following example. Consider the connectivity graph given by the couplings in matrix A to be the one shown in Figure 3.9, where each box represents a different subdomain. Blue nodes are those belonging to the first slice of their linelet. In the first step, local slices are assembled. This is shown between figures 3.9a and 3.9b and was previously shown in Figure 3.1. Also, red dotted lines in Figure 3.9 denote the couplings between different subdomains (i.e., couplings between local slices and the halo). The process of coupling different slices is performed iteratively via successive halo updates. Starting from those slices candidates to be the first of their linelet (depicted as slices containing blue nodes in Figure 3.9, these blue nodes may represent, for example, the nodes or elements adjacent to the airfoil in Figure 1.7), if any of these slices are coupled to the halo, each one of them is joined with the corresponding node in their respective subdomain. This process is iterated Np− 1times until at most Np slices are joined together or until there are no more slices to join. Figures 3.9c to 3.9f show how the linelets grow to cross one interface per step. Hence, the iterative process shown in Figure 3.9 can be summarized as follows: 1. Each subdomain builds its own set of slices as described in section 3.1 (and further developed in 3.1). It is also important that the information about their couplings to the halo (red dotted lines in Figure 3.9) is stored. 2. To grow the linelets, we must select the slices that will be the first of their linelet. In this example, they are taken as such those including blue nodes in Figures 3.9a to 3.9f. As previously mentioned, they may refer to the nodes or elements adjacent to the airfoil in Figure 1.7. 67 3.3. Growing the linelets: joining slices from different subdomains (a) Sketch of a hybrid mesh divided into several subdomains. Blue nodes indicate that their slice is a candidate to be the first in their respective linelet. (b) Each subdomain constructs its slices independently from one another. When a slice has a coupling with the halo, this coupling is stored to be used when growing the linelets. Red dotted lines denote couplings with the halo. (c) After a first halo update, linelets are grown by connecting the first slices (those containing blue nodes) with the slices they are coupled to through the halo. (d) A second halo update is performed to incorporate three slices into the linelets by joining the previously grown linelets with the slices they are coupled to via the halo. (e) A third halo update is used to grow the linelets to include four slices. (f) After the final halo update, any coupling not included in a linelet is discarded. Figure 3.10 shows this final stage. Figure 3.9: Example of the linelet construction process using slices from different subdomains. 68 3.3. Growing the linelets: joining slices from different subdomains Figure 3.10: Final stage after assembling the slices and growing the linelets. Only couplings within the linelets are kept, and every other coupling is discarded. 3. An halo update operation is repeated Np− 1times ( Np = 5 in the previous example). By this operation, each slice knows to which other slice is coupled (change from dotted to continuous lines in Figure 3.9c). In this step, the communication scheme to be used in section 2.2.2 is built. 4. After Np− 1iterations, if any remaining coupling has not been used to join slices, that coupling is discarded. Also, if no linelet has grown from any given boundary node, this node is treated as a node not included in any linelet. To apply the Schur Complement method to solve equation (1.39) for each node included in a linelet, we need to identify the interface nodes, which we defined to be the first node of each local linelet after the interface (by interface, it is meant the boundary between different subdomains). Figure 3.10 shows these as empty dots. This way, every inner node is coupled to other inner or interface nodes, just as the Schur algorithm requires. Also, note that the first slice of the linelet does not have interface nodes. These nodes are taken to be the first one of the i−th linelet slice for i > 1. After this step, each process knows to which linelet its slices belong and their position in their respective linelet. The next step is to define the communication scheme. To solve the linear system for the interface nodes of each linelet (equation (1.47) ), we assign each linelet to a specific process. This process solves equation (1.47) for the corresponding linelet. Hence, the first step is to map each linelet to a given process. In our case, we opted to distribute the linelets as evenly as possible. Also, within this approach, linelet loops are naturally avoided. This is the method for growing the linelets we have developed and published in the work published in [2]. However, the linelet growth method has been extended and generalized, not requiring a starting point specification. This new approach is explained in the next section and the one we will continue to use in 69 3.3. Growing the linelets: joining slices from different subdomains the present work. For a more detailed description of this approach, the reader can refer to Appendix B. 3.3.2 Growing the linelets without specifying a starting point An alternative approach adopted in this work involves growing linelets without specifying a starting slice. This method is more general and removes the need for user input. However, attention must be taken to prevent the formation of linelet loops. Consider the example in Figure 3.11, where a slice with global ID j is coupled on one side to a slice with global ID i and on the other to a slice with global ID k . During a halo update communication step, coupling information is exchanged. Then, the process managing slice i has the information about its coupling with slices jand k, and the same with the process managing slices k and j. Subsequent halo updates further propagate this information. For instance, the process handling slice j informs the process managing slice k of its coupling with slice i . As a result, the process responsible for slice k now has the information of both its direct coupling with slice j and the indirect coupling between slices j and i and similarly with the process managing slice i . At this point, processes containing i , j , and k have the complete information about the global structure of the linelet they belong to. Because the linelet IDs are globally unique, the coupling information also includes the ranks of the processes managing the coupled slices. This guarantees that all data necessary for constructing the communication scheme is readily available. This process is illustrated in Figure 3.11. Therefore, this approach begins by assigning global IDs to all slices coupled to the halo. To grow the slices, consider the example in Figure 3.12, where ten slices are coupled and numbered as shown. The first step involves performing an allgather communication across all processes to share information about the number of slices each process has. This allows every process to assign a unique global ID to its slices. To do this, each process assigns a variable called nSlices, corresponding to its number of slices. Then, an allgather communication is performed with all the processes. This way, each process knows the number of slices of every other process. Then, each process assigns a global ID to each slice, ensuring no ID is repeated. For example, if the array resulting from this allgather communication is [3, 2, 4, 1, 5], then the slices of the first process will be assigned global IDs 1, 2, and 3, the slices of the second process will be assigned global IDs 3 and 4, and so on. Note that there is no need to assign an ID to the slices that are not coupled to the halo. Within our implementation, we will assign a global ID to every slice, even those not coupled to the halo, but with a caveat: slices not coupled to the halo will have a global ID equal to -1, while positive slices are reserved only to those coupled to the halo. Also, 70 3.3. Growing the linelets: joining slices from different subdomains Figure 3.11: Example of the communication scheme used to grow the linelets without specifying a starting point. slices with only one coupling to the halo will have the lowest ID. Later, we will see that this will help prioritize linelets growing from the boundary layer. This global numbering of the linelets is done by doing the aforementioned allgather but, instead of nSlices, with another variable nSlicesCoupled, which stores the number of slices coupled to the halo. By the end of this step, every process containing a slice will have an array, which we called IdGlobal of size nSlices containing the global ID of each slice: local slice i has a global ID of IdGlobal[i]. Note from Figure 3.11 and the explanation provided that the information flows in two directions, and therefore, two different ‘channels of communication’ are used. The main idea is to assign to each slice a two-dimensional array, with each dimension corresponding to a communication channel through which information flows in a specific direction. This is the basis of this process. Consider the example in Figure 3.12, which depicts a set of ten coupled and numbered slices. The upper section of the Figure illustrates the coupled nodes forming slices and their corresponding global IDs, while the lower section displays the entries of the AuxField array assigned to each slice. The halo entries, to which the nodes at the extreme of the slices are coupled, are represented as dotted entries. In this example, the first degree of freedom facilitates the flow of information from left to right, while the second handles the flow from right to left. Therefore, the process begins by defining an auxiliary field with two degrees of freedom, called AuxField, where each node is mapped to two entries: node i is mapped to AuxField [ i, 1] and AuxField [ i, 2]. This array serves as the basis for successive halo updates. For every node located at the extremity of a slice, the corresponding entries in AuxField are filled with the global ID of their slice. Although all entries associated with the slice could be filled, only the entries required for communication during the halo update are necessary, as only they influence the result. We must now define a parameter Np , which represents the number of times 71 3.3. Growing the linelets: joining slices from different subdomains Figure 3.12: Example of the variables used to grow the linelets without specifying a starting point. this process will be iterated. In the following example, we consider Np = 2, as this choice highlights specific issues that must be addressed with this approach. Figure 3.13 illustrates the flow of information over two halo communication steps for the case depicted in Figure 3.12. The colored entries in the AuxField variable indicate relevant entries, while all others are byproducts of the method and are not essential. In the first step, the entries in the AuxField array associated with the nodes at the extremes of the slices are filled with their corresponding slice’s global ID. Then, a halo update is performed. The subsequent step is to update the entries of the AuxField array associated with the nodes at the extremes of the slices. In this step, it is important to be consistent with the direction in which the information must flow. For example, consider the slice with global ID 7 (the second slice from the left in Figure 3.13). After the first halo update, if the left node receives an input with a value of 4, it must be propagated to the right node to pass it to the adjacent slice (to the right) in the next halo update. Similarly, if the right node receives a value of 9, this value must be copied to the left node for propagation to the adjacent slice (to the left) in the subsequent halo update. For a better understanding of the process, in the diagram of figure 3.13, we have kept that the flow of information to the right goes through the first channel while the flow to the left goes through the second one. However, this is not necessary. It is sufficient that, for example, the right node is updated with both inputs corresponding to the two degrees of freedom of the left halo node. If it receives two entries of different values, the relevant information is in the one whose value it has not yet received. Special attention is required when dealing with single-node slices. In such cases, the relevant input (i.e., the one not previously received) must be copied from each halo node into each degree of freedom of the inner node. While this step might initially seem confusing, it will become clearer with the explanation of the next step. After this step, each slice has information about its couplings with adjacent slices and their respective IDs. For example, slice 4 ‘knows’ is coupled with slice 7 on the right, slice 7 ‘knows’ is coupled with slice 4 on the left, slice 9 on the right, and so on. At this point, we introduce a new variable that will be useful moving forward: an array called CoupledSlices, with dimensions 2( Np + 1) ×nSlices , where nSlices is the number of slices in the current process, and the array has 72 3.3. Growing the linelets: joining slices from different subdomains Figure 3.13: Detailed example of the process used to grow the linelets without specifying a starting point. Figure 3.14: Detailed example of the process used to grow the linelets without specifying a starting point. only one degree of freedom. Each slice in this array corresponds to 2( Np + 1) entries in CoupledSlices. The entries corresponding to the i -th slice of the process are located in CoupledSlices[2( Np +1)(i-1)+1:2( Np +1)i]. We assign the global ID of the slice to the first entry: CoupledSlices[2( Np +1)(i-1)+1] = IdGlobal[i], and initialize the remaining entries to zero. This array will store the order in which the slices are coupled. Figure 3.14 illustrates the evolution of the CoupledSlices section corresponding to each slice for the example shown in Figure 3.12. After the first halo update, CoupledSlices is updated according to the couplings received. For example, since slice 4 received only the coupling with slice 7, it adds this value to CoupledSlices. Similarly, slice 7 received the coupling with slice 4 and slice 9, so these two values are added. The coupling received from the left is placed to the left of 7 (at the first entry, i.e., at the 2( Np +1)(i-1)+1 entry), and the coupling from the right is added to the right 73 3.3. Growing the linelets: joining slices from different subdomains Figure 3.15: Example of the evolution of the CoupledSlices array from the example shown in figure 3.12 after the first halo update. of 7 (i.e., replacing the first zero in CoupledSlices[2( Np +1)(i-1)+1:2( Np +1)i]) in CoupledSlices. It is important to note that although it does not matter whether couplings from the left or right are placed on the left or right side of CoupledSlices, consistency must be maintained, and the same approach should be followed after each halo update. The idea is that any new coupling information is added to CoupledSlices. However, if a coupling value already exists in CoupledSlices, it is not added again, ensuring that no entries are repeated for a given slice. The third column from the right in Figure 3.14 shows the state of CoupledSlices after the second halo update is completed. To build the communication scheme correctly, each slice must be aware of its position in the linelet to which it belongs, which leads us to the next step. Now that each slice knows its couplings and those of its neighboring slices, the next step is to identify the different linelets and determine each slice’s position within its respective linelet. It is important to note that the CoupledSlices array, as shown in the third column of Figure 3.14, does not yet reflect the global structure of the linelets. For example, we might initially think that slice 4 (in the first row) belongs to the linelet formed by slices 4, 7, and 9, while slice 7 belongs to the linelet formed by slices 4, 7, 9, and 2. This would create an inconsistency in the linelet structure, potentially disrupting the solution of the corresponding linear system. To tackle this, an additional step is needed to ensure the consistency of the information associated with each process slice. In this step, each slice identifies the minimum value in its respective CoupledSlices section, i.e., the slice with the smallest global ID to which it is coupled. This value is then communicated to the neighboring slices through a halo update, similar to the previous process involving AuxField. Note that only one degree of freedom is required for this step. An example of this process is illustrated in Figure 3.12. Now, suppose after this halo update, it does not receive from its neighbor 74 CHAPTER 4 Implementation and performance Part of the contents of this chapter have been published as: R. de Olazábal, R. Borrell and O. Lehmkuhl. ‘An algebraic global linelet preconditioner for incompressible flow solvers’. In: Journal of Computational Physics 514 (2024), and disseminated in the following conferences: the 15th JLESC Workshop, Bordeaux, France, 2023, and the 1st Math 2 Product (M2P) conference for Emerging Technologies in Computational Science for Industry, Sustainability and Innovation, Taormina, Sicily, 2023. In this section, we present numerical results and compare the performance of the diagonal preconditioner, the LLP, i.e., the linelet preconditioner when communications are not allowed, and the Global Linelet Preconditioner (GLP, i.e., when communications are allowed). It has been seen in previous works [6] that for anisotropic meshes with highly stretched elements, LLP reduces the number of PCG iterations. Also, we observed that when this anisotropy is due to mesh stretching, GLP significantly reduces the number of iterations in the solution of the pressure equation and, consequently, the total CPU time. To evaluate the proposed method, we solved Poisson’s equation using the finite element method in a large-scale, unsteady, incompressible flow problem with a Reynolds number of Re = 10 5 . The Reynolds number or the choice of turbulence modeling approach affects only the r.h.s.. In our case, we have obtained the r.h.s. for LES. However, we expect a very similar behavior if we had used RANS. It is important to highlight that RANS can be used in meshes with higher anisotropy than those usually employed in LES. This is particularly advantageous for the proposed method, as it has demonstrated faster convergence in cases with higher aspect ratios. It is important to clarify that we focused on solving Poisson’s equation for the pressure correction rather than the complete set of governing equations in addressing these challenges. As such, the ensuing results correspond exclusively 81 4.1. Implementation to the solution of the pressure correction equation. 4.1 Implementation The present solver has been developed within Alya, the multi-physics finite element simulation software designed and maintained by the Barcelona Supercomputing Center (BSC). Alya is optimized for high-performance computing (HPC) environments, ensuring scalability on both CPUs and GPUs through hybrid parallelization models that combine MPI, OpenMP, and OpenACC. The software has undergone rigorous verification, validation, and optimization, demonstrating its reliability in solving complex fluid dynamics problems, among other physical phenomena [7,1,3,4,5]. Furthermore, Alya is one of the twelve simulation codes included in the Unified European Applications Benchmarks Suite (UEABS), underscoring its adherence to the highest HPC standards. These attributes make it a suitable platform for implementing and testing the present solver in scientific and engineering applications. Alya follows a modular architecture consisting of a core kernel and multiple specialized modules. The kernel provides fundamental functionalities for solving discretized PDE, while the modules define the specific physics of a given problem. The software supports a diverse range of physical simulations, including incompressible and compressible CFD, computational solid mechanics (CSM), heat transfer through diffusion, convection, and radiation, particle transport, electrophysiology, reactive species transport, and nuclear reactions. Depending on the multi-physical nature of a given problem, these modules can be coupled to enable complex simulations, such as spray combustion modeling, which integrates turbulent flow, reacting chemical species, heat transfer, and evaporating droplets. 4.2 Workflow The workflow of Alya follows a standard structure used in engineering solvers for nonlinear PDEs. The process begins with a preprocessing stage in which the input data is read, memory is allocated, communication patterns are established, and data structures are initialized. Next, a time-stepping loop is used for transient problems. For nonlinear systems, an additional iterative loop runs at each time step to assemble and solve the nonlinear algebraic system, repeating until the error falls below a set threshold. Once the error falls below a specified threshold, the simulation moves to the next time step. After all time steps have been completed, the post-processing phase deallocates memory and ends the simulation. 4.3 Computational Resources All simulations presented in this thesis were executed on the MareNostrum V supercomputer at BSC. As a pre-exascale EuroHPC system, MareNostrum V offers a peak computational power of 314 petaflops, leveraging a combination of Bull Sequana XH3000 and Lenovo ThinkSystem architectures. The system 82 4.4. Numerical verification of the solver is divided into four specialized partitions to cater to various HPC workloads. The primary partitions include: • The General Purpose Partition (GPP), delivering 45 petaflops with 6480 standard nodes and 72 high-bandwidth memory (HBM) nodes, powered by Intel Sapphire Rapids processors. • The Accelerated Partition (ACC), offering 230 petaflops through 1120 nodes equipped with Nvidia Hopper GPUs and Intel Sapphire Rapids processors. Additional partitions featuring Nvidia Grace CPUs and next-generation accelerators will further extend the system’s capabilities. MareNostrum V employs a high-performance fat-tree network topology and an advanced storage system providing 248 petabytes of capacity, with read/write speeds of 1.6 TB/s and 1.2 TB/s, respectively. A long-term archive supplements this with an additional 402 petabytes of tape-based storage. 4.4 Numerical verification of the solver To ensure accuracy, reliability, and robustness, the solver underwent a comprehensive verification and validation process using multiple test cases. Its performance and scalability were also tested across different mesh sizes and element types. In this section, we present numerical results and analyze the performance of three preconditioners: the diagonal preconditioner, the LLP, which operates without inter-domain communication, and the GLP, which allows for communication between subdomains. Previous studies [6,2] have shown that for highly anisotropic meshes with highly stretched elements, LLP effectively reduces the number of PCG iterations. Moreover, when the anisotropy results from mesh stretching, GLP has been observed to further accelerate convergence in solving the pressure equation, significantly reducing total computational time. To evaluate the proposed approach, we applied it to solve Poisson’s equation using the finite element method within a large-scale, unsteady, incompressible flow simulation at a Reynolds number of Re = 10 5 . Notably, when solving the pressure equation, the resulting matrix corresponds to a Laplacian operator that depends solely on the mesh structure rather than the flow properties. The Reynolds number and the turbulence modeling approach influence only the r.h.s. of the equation. In our case, this term was obtained from a LES. However, we anticipate that similar behavior would be observed if a RANS model were used. It is worth emphasizing that RANS is typically applied to meshes with even higher anisotropy than those used in LES, making it particularly well-suited for the proposed method, as its convergence benefits become more pronounced in cases with greater aspect ratios. Furthermore, we tested the solver using a random vector on the r.h.s. and observed a similar convergence behavior. 83 4.4. Numerical verification of the solver It is important to clarify that we focused on solving Poisson’s equation for the pressure correction rather than the complete set of governing equations. As such, the following results correspond exclusively to the solution of the pressure correction equation. It is also worth mentioning that, in our experience, the linelet preconditioner, designed to exploit the mesh anisotropy in its favor, has demonstrated improved performance with increasing aspect ratios. Specifically, a higher aspect ratio correlated with a reduced number of iterations required for convergence, indicating faster convergence rates (we define the aspect ratio as the geometrical aspect ratio, representing the ratio of the maximum to minimum length of a given element). Consequently, the presented results are particularly advantageous for scenarios characterized by higher aspect ratios. The simulations for this case were executed utilizing 48 to 768 processes, with the parameter Np (previously discussed in chapter 3) selected for each case to ensure the inclusion of every linelet slice. Importantly, as the number of processes increased, all problem parameters remained constant, allowing us to conduct a robust, strong scalability test. In order to study the convergence of the GLP, we compare the value of the normalized residual given by Res =||A.xi−b|| ||b|| ,(4.1) for GLP, LLP , and the diagonal preconditioner, where xi is the value of the increment in pressure at the i−th iteration. Since we are solving for the increment in pressure, we take the initial value of the unknown (x0) to be  0. In the following analyses, we employ three different meshes to study the proposed generalization of the linelet preconditioner. These include a hybrid mesh comprising 1.64 million elements, with 60% of them located within the boundary layer and a maximum aspect ratio of 17 . 05. Additionally, we utilize a structured mesh with 0.98 million elements, all confined within the boundary layer and with the same maximum aspect ratio of 17 . 05. To explore the efficacy of the proposed improvements under more challenging conditions, we introduce an extreme case of a structured mesh with 5 million elements. All elements are located within the boundary layer in this configuration, and the maximum aspect ratio reaches 284.68. Figure 4.1 shows the convergence of the residual throughout each iteration of the Preconditioned Conjugate Gradient (PCG) method for two of the analyzed cases: 1.64 and 0.98 million element meshes. Note how in Figure 4.1a there is no significant difference in the convergence during the first steps, but there is a threshold after which GLP accelerates its convergence. On the other hand, when all the elements are located within the boundary layer, which is the case depicted in Figure 4.1b, it can be seen that the improvement in convergence is almost immediate. If we instead plot the number of iterations required for each method to reduce the residual by six orders of magnitude (i.e., a residual such that the 84 4.4. Numerical verification of the solver (a) Case run on a mesh with 1.64 million elements, with 60% in the boundary layer. GLP reaches a residual of 10 −6 1.62 times faster than LLP and 4.54 times faster than the diagonal preconditioner. (b) Case run on a mesh with 0.98 million elements entirely located within the boundary layer. GLP reaches a residual of 10 −6 2.97 times faster than LLP and 9.27 times faster than the diagonal preconditioner. Figure 4.1: Comparison of PCG convergence using GLP, LLP, and diagonal preconditioners in PCG Method on meshes with varying element counts and boundary layer coverage. Residual is plotted against the number of iterations. Simulations executed on 384 processes. 85 4.4. Numerical verification of the solver quotient between the residual at the i-th iteration and the RHS equals 10 −6 , see equation (4.1)), we obtain the graphs shown in Figure 4.2. It is important to note that GLP requires only 534 iterations to reach the desired tolerance, whereas the LLP requires at least 682 iterations, and the diagonal preconditioner requires between 3043 and 3046 iterations. Therefore, the proposed method reduces the number of iterations by a factor of at least 1.27 compared to the LLP and a factor of 5.69 compared to the diagonal preconditioner. It is also noteworthy to mention that both LLP and GLP exhibit notably accelerated convergence compared to the scenario where only 60% of elements are within the boundary layer (see Fig. 4.2a). On the other hand, the convergence for the LLP worsens when the number of processes increases. Furthermore, this improvement has been observed to increase for higher aspect ratios or meshes with a larger number of elements. However, even though the GLP takes a lower amount of iterations to converge to a desired tolerance, it should also be taken into account that each iteration takes more CPU time since more operations (mainly communications) are required. Thus, for the cases previously shown in Figure 4.1, Figure 4.3 shows the convergence of the residual at every PCG iteration vs. execution time. Figure 4.3a shows the convergence in time for the case previously shown in Figure 4.1a, and Figure 4.3b shows the same relation for the case previously shown in Figure 4.1b. Note how in these cases the GLP curve aligns closely with LLP due to the longer duration of each GLP iteration. Additionally, the data presented in Figures 4.3a and 4.3b can be expressed in terms of execution time required to reach a specific tolerance. This is demonstrated in Figures 4.4a and 4.4b, which show the time per iteration for each preconditioner. Figure 4.4a indicates that the diagonal preconditioner has the shortest time per iteration, followed by LLP, and finally, GLP has the longest time. However, the total time required to achieve a certain tolerance is determined by the time per iteration and the number of iterations needed. In most cases, the number of iterations required follows an inverse relation, with GLP requiring the least number of iterations, followed by LLP, and finally, the diagonal preconditioner. As a result, these two factors determine which is the fastest preconditioner. Figure 4.4b shows the total time required by the PCG solver to converge to a residual of 10 −6 . Note that GLP takes less time to achieve the desired solution since iterations are significantly reduced in this case. It is important to note that, for all the studied cases, GLP always reduced the number of iterations required to reach a certain tolerance. Nevertheless, this does not always mean an improvement in time. Let us show other cases where GLP represents little to no improvement in time and one where it does represent a significant improvement. Figure 4.5a compares the residual against the number of iterations for the three preconditioners under study: Diagonal, LLP, and GLP, for a hybrid mesh of 1.64 million elements where 40% of them are within the boundary layer and using 384 processes. Note how GLP improves the convergence: For all the studied cases, it always takes the lowest number of iterations. Nevertheless, if we compare the time required to converge, we see that there is little to no improvement at all. This is shown in Figure 4.5b. Hence, even though GLP converges in a lower amount of iterations, 86 4.4. Numerical verification of the solver (a) Case run on a mesh with 1.64M elements, 60% of them within the boundary layer. (b) Case run on a mesh with 0.98M elements, all within the boundary layer. Figure 4.2: Comparison of the number of iterations required to reach a tolerance of 10 −6 for GLP, LLP, and diagonal preconditioner on meshes with varied element counts and boundary layer coverage. 87 4.4. Numerical verification of the solver (a) Convergence of the residual in time for the case shown in Figure 4.2a. GLP achieves a residual of 10 −6 1.17 times faster than LLP and 2.07 times faster than the Diagonal preconditioner. (b) Convergence of the residual in time for the case shown in Figure 4.2b. GLP achieves a residual of 10 −6 2.12 times faster than LLP and 4.08 times faster than the Diagonal preconditioner. Figure 4.3: Comparison of the residual convergence in time between GLP, LLP, and diagonal preconditioner. Simulations executed on 384 processes. 88 4.4. Numerical verification of the solver (a) The diagonal preconditioner exhibits the shortest time per iteration, followed by LLP. As expected, GLP takes the longest time per iteration. (b) Time required to reach a tolerance of 10 −6 for the GLP, LLP , and diagonal preconditioner. GLP requires less time due to its faster convergence. Figure 4.4: Comparison of iteration time and total convergence time for GLP, LLP, and diagonal preconditioners in achieving 10−6tolerance. 89 4.4. Numerical verification of the solver (a) Residual vs. the number of iterations. GLP achieves convergence in the lowest number of iterations among all tested preconditioners: 10 −6 1.40 times faster (in terms of the number of iterations) than LLP and 3.80 times faster than the Diagonal preconditioner. (b) Convergence of the residual over time. In this case, GLP offers little to no improvement over the LLP. GLP reaches a residual of 10−61.07 times faster than LLP and 1.81 times faster than the Diagonal preconditioner. Figure 4.5: Residual convergence comparison in PCG iterations and time: GLP, LLP, and diagonal preconditioner for a mesh of with 1.64 million elements of 40% of them located within the boundary layer. The simulation was run in 384 parallel processes. the time required to converge is almost the same as in the LLP case. These observations emphasize the trade-off between iterations and time for different preconditioning methods. To further understand the behavior of the proposed method, if we run the same case as before but this time in 768 processes, we have the results shown in Figure 4.6. Figure 4.6a shows the same behavior as the one shown in Figure 90 4.4. Numerical verification of the solver (a) (b) Figure 4.10: Improvement in the preconditioned residual convergence time for different tolerances and percentages of elements within the boundary layer. The color scheme is the same as in Figure 4.9 97 4.5. Analysis of performance with the number of partitions, degrading the preconditioner and consequently increasing the number of iterations required to reach a certain tolerance. This is consistent with our experience, where we have consistently observed that GLP is the preferred option for high anisotropic meshes with a high percentage of elements within the boundary layer. Two main considerations should be noted: • Based on our experience, a "high" percentage of elements within the boundary layer is typically considered to be 40%. The time improvement becomes more significant for meshes with more than 40% of elements in the boundary layer. Beyond 60%, the improvement is even more pronounced. If the expected improvement is small, implementing GLP may not be worthwhile. For example, in Figure b) of Figure 4.9, for a mesh with 50% of the elements in the boundary layer and a tolerance of 10 −3 , GLP requires 92.1% of the time needed by LLP. In this case, the improvement might not justify the additional computational effort. • The number of partitions is not the only factor influencing performance— the quality of the mesh partitioning relative to the linelets also plays a crucial role. For example, in a high-partition scenario, the mesh may be partitioned such that LLP discards fewer couplings, resulting in faster convergence. Conversely, for low-partition scenarios, if the mesh is partitioned so that the number of couplings that LLP discards is high, this can lead to slower convergence. These observations support the idea that the proposed method presents a promising alternative to the state-of-the-art linelet preconditioners, particularly for high anisotropic meshes with a large proportion of elements in the boundary layer. However, the improvement in time is most notable for high tolerances and meshes with a significant boundary layer percentage. Additionally, the number of partitions and the quality of mesh partitioning are key factors that influence the solver’s performance. 4.5 Analysis of performance This section presents an analysis of the performance of the preconditioner. We start by analyzing the execution of the solver by inspecting the traces. This provides a clear representation of the workload distribution and the overhead introduced by the proposed method. We then study the performance of the proposed method by using standard performance metrics. The analysis focuses on speed-up, load balance, communication efficiency, and overall parallel performance. 4.5.1 Execution Analysis The execution and implementation of the proposed method were analyzed using Extrae, a performance analysis tool used to probe and register parallel executions. Extrae generates trace files that can be visualized using Paraver, a 98 4.5. Analysis of performance Figure 4.11: Trace of the parallel execution, showing the time spent in each step of the PCG for the second iteration of the first time step. performance analysis and visualization tool. Both tools were developed at the Barcelona Supercomputing Center. The analysis was conducted on a cavity flow problem with a mesh of 340000 elements, with 50% of them within the boundary layer. The simulation was executed using 32 processes, with a parameter Np = 100, such that every linelet slice is included in the global linelet structure. The slices were grown with the algebraic method and a parameter α = 2. The results for the second iteration of the first time step are shown in Figure 4.11. For the PCG operations, we use the same notation as in section 1.5.3. A key aspect to note is that the time spent in SpMv and dot product computations includes the communication required for their execution. This explains why the dot product rTs takes significantly longer for certain processes: while some have already completed the preconditioning step, others are still in the process of finishing it. Step 6 (communication of interface values) is subdivided into two phases: mapping data to the sending buffer and performing the actual communication. Additionally, in Alya, the first thread corresponds to the master and does not participate in preconditioner computations. Therefore, it is excluded from the present analysis. 99 4.5. Analysis of performance Figure 4.12: Zoomed-in trace of the preconditioner execution for the second iteration, corresponding to the white rectangle in Figure 4.11. Another important observation is the load imbalance between processes 2–16 and 17–32, which results from the non-uniform distribution of the boundary layer. Processes handling a larger portion of the boundary layer contain a higher concentration of linelets. If a process has linelets, it applies the linelet preconditioner; otherwise, it resorts to the diagonal preconditioner. Since the linelet preconditioner requires more computations than the diagonal one, processes with a higher density of linelets take longer to complete the preconditioning step. To further illustrate this, Figure 4.12 provides a zoomed-in view of the white rectangle in Figure 4.11, highlighting the time spent in each step of the preconditioner. The figure contrasts the execution of a process with linelets (thread 16) against one without linelets (thread 17). From Figures 4.11 and 4.12, it can be seen that the workload associated with solving for the interface (step 5) is evenly distributed across all processes. However, after the communication step (step 6), thread 17 proceeds directly to the diagonal preconditioner, while thread 16 solves the linear system associated with the inner nodes of its slices. In our approach, this last step is done by iterating over the slices. During each iteration, the process computes the r.h.s. bp−Apsxs (step 7), solves for the inner nodes (step 8), and maps the solution back to its original ordering (step 10). Figure 4.12 also provides a more detailed view of the iterative process used to solve for the inner nodes within the slices. 100 4.5. Analysis of performance 4.5.2 Scalability Analysis and Performance Metrics Given the complexity of the problems addressed in this work, computational efficiency is critical to minimizing energy consumption, costs, and processing time. Thus, we must ensure that the implementation scales well with the number of CPU cores allocated for the task. To this end, we analyze different performance metrics of the solver. Let Npref denote the minimum number of cores used in a performance test, with the corresponding elapsed time given by Tref . If the same problem is executed on Np cores with an elapsed time of TNp , the relative speed-up is defined as: SU(T) = Tref TNp (4.2) In our case, TNp is considered to be the maximum elapsed time across all processes. We also consider the speed-up of the mean execution time across all processes, which is given by replacing TNp by TNp = ⟨Ti⟩Np i=1 in equation (4.2) : SU(⟨T⟩) = Tref ⟨Ti⟩Np i=1 (4.3) However, a high speed-up does not guarantee efficient resource utilization. To assess this, three key metrics provide a more in-depth analysis: load balance, communication efficiency, and parallel efficiency. Load balance quantifies how evenly the computational load is distributed among different processing units, therefore reducing bottlenecks due to idle processes. It is computed as the ratio between average useful computation time across all processes, Tw = ⟨Ti w⟩Np i=1 , and the maximum useful computation time (also across all processes), Tmax w= maxNp i=1 Ti w: LB =Tw Tmax w (4.4) On the other hand, communication efficiency is the maximum, across all processes, of the ratio between useful computation time and total runtime. It measures the proportion of time spent on computations relative to communication overhead. It is defined as: CE =Tmax w TNp (4.5) Poor communication efficiency identifies cases where excessive time is spent on communication rather than executing useful computations. Two primary factors influence this metric: 1) synchronization delays, where processes remain idle at communication points due to workload imbalance, and 2) data transfer overhead, which arises when processes exchange large volumes of data relative to the network’s network capacity or with high-frequent communication. Finally, parallel efficiency assesses the overall utilization of computational resources to achieve this speed-up. It is expressed as: 101 4.5. Analysis of performance PE =Tw TNp =LB ×CE (4.6) Note that the parallel efficiency can be computed as the product of load balance and communication efficiency. Thus, the importance of minimizing communication overhead and ensuring a balanced workload to utilize resources in parallel computing efficiently. In order to evaluate the performance of the solver, we conducted a series of tests on a large-scale, unsteady, incompressible flow problem with a Reynolds number of Re = 10 5 . The simulations were executed using 56 to 1344 processes, with the parameter Np (discussed in section 3) selected for each case to ensure the inclusion of every linelet slice. Importantly, as the number of processes increased, all problem parameters remained constant, allowing us to conduct a robust, strong scalability test. For the following results, we used a cavity flow problem with meshes of 2.25 million elements with 40% to 90% of elements within the boundary layer. Note that more elements within the boundary layer imply 1) larger linelets and/or 2) more linelets. Since the main aspect of the solver we are discussing is the implementation of the newly developed preconditioner, in this section, we only consider computation and communication times spent within the preconditioner. Figures 4.13 and 4.14 show the parallel performance analysis for the cavity flow problem with 40% of elements within the boundary layer. The results are presented as a function of the number of processes, with the following metrics: speed-up, meantime speed-up-, load balance, and parallel efficiency. Ideally, the speed-up should align with the reference dashed red line, which represents perfect linear scalability. Any deviation from this ideal behavior indicates inefficiencies stemming from communication overhead and workload imbalance. Similarly, communication efficiency should ideally remain close to 1, but it naturally degrades as inter-process communication costs increase. A perfectly balanced workload would yield a load balance value of 1, while lower values indicate disparities in computational workload distribution across processes. Finally, parallel efficiency integrates the effects of speed-up, communication efficiency, and load balance, offering a comprehensive measure of how effectively the preconditioner scales across multiple processes. Table 4.1 provides a comparison of the LLP and GLP performance metrics. The results show that for most cases, the GLP outperforms the LLP in terms of speed-up, meantime speed-up, and load balance while maintaining a similar parallel efficiency. The results for the diagonal preconditioner are shown as a reference for comparison. The main reason for the speed-up improvement in the GLP over the LLP is that the work is spread more evenly among processes, as shown by the load balance metric. This is particularly important in the context of the present solver, as the preconditioner is usually the most computationally expensive part of the simulation. In the LLP, the processes containing slices must solve a TDMA, while the rest of the processes apply a dot product related to the diagonal preconditioner. This leads to idle processes and a poor load balance. In contrast, the GLP allows all processes to distribute the work more evenly 102 4.5. Analysis of performance since not only the processes containing slices need to solve a TDMA, but also the processes in charge of solving for the interface nodes need to apply a TDMA. However, for a large number of processes, the communication overhead may become a bottleneck, as evidenced by the decrease in the speed-up and parallel efficiency with respect to the LLP case. Also, note that the quality of the mesh partition is crucial for the performance of the GLP, as the communication scheme is determined by the linelet partitions, which are determined by the mesh partition. Therefore, a poor load balance in the linelet preconditioner (LLP or GLP) can be a consequence of a mesh partition that cuts the linelets in a way that the slices are not evenly distributed among the processes. This is the reason why we compare the LLP and GLP performance since the GLP is a generalization of the LLP that allows communication between slices. Therefore, the fact that the GLP outperforms the LLP in most cases is an indicator of 1) a good implementation of the GLP and 2) that the present solver is able to take advantage of the GLP communication scheme. Finally, note that since the diagonal and the LLP does not require communications, their communication efficiency is always 1, and therefore their load balance and parallel efficiency are the same. Pr LLP GLP SU(T)SU(⟨T⟩)LB/PE SU(T)SU(⟨T⟩)LB CE PE 56 1.00 1.00 0.42 1.00 1.00 0.44 0.89 0.95 112 1.55 1.82 0.35 1.80 1.92 0.41 0.90 0.91 168 2.31 2.60 0.37 2.52 2.97 0.38 0.91 0.89 224 2.70 3.29 0.35 3.68 3.92 0.40 0.94 0.84 336 3.56 4.60 0.30 4.09 5.89 0.33 0.90 0.81 448 4.02 5.65 0.30 5.07 7.85 0.33 0.88 0.74 672 5.84 8.29 0.29 6.75 11.90 0.32 0.87 0.76 896 7.20 10.26 0.29 7.62 15.92 0.30 0.86 0.64 1344 9.20 14.43 0.27 9.24 23.85 0.29 0.78 0.52 Table 4.1: Comparison of LLP and GLP performance metrics for a cavity flow problem with 2.25 million elements, 40% of them within the boundary layer. Pr refers to the number of processes. A more detailed analysis of the GLP ’s performance is presented in Figures 4.15, 4.16 and 4.17, where we also include the communication efficiency. It is observed that communication efficiency decreases as the number of processes increases, with this effect being more pronounced for meshes containing a higher percentage of nodes within the boundary layer. This behavior is attributed to the rising communication overhead, which escalates with the number of processes and the number and length of the linelets. As the number of interface nodes to be communicated grows, so does the associated overhead. This is also reflected in the minimal to no improvement in speed-up and the poor parallel efficiency observed when running on a large number of processes. However, it is important to highlight that the load balance remains con103 4.5. Analysis of performance (a) Speed-up (b) Mean time speed-up Figure 4.13: Speed-up as a function of the number of processes for a cavity flow problem with 2.25 million elements, 40% of them within the boundary layer. 104 4.5. Analysis of performance (a) Load balance (b) Parallel efficency Figure 4.14: Load balance and parallel performance as a function of the number of processes for a cavity flow problem with 2.25 million elements, 40% of them within the boundary layer. 105 4.6. Preprocessing stage (a) Speed-up Figure 4.15: Speed-Up for the GLP as a function of the number of processes for a cavity flow problem with 2.25 million elements Different percentages of elements within the boundary layer. sistently maintained across all numbers of processes, improving with a higher percentage of elements within the boundary layer. Two key factors explain this: 1) as the percentage of elements within the boundary layer increases, so does the number of nodes within the linelets, and 2) the communication strategy—where each linelet is assigned to a specific process to handle its interface nodes—results in a more even distribution of the load as the number of linelet nodes increases. 4.6 Preprocessing stage The preprocessing stage is the first step in the solver’s workflow. During this stage, the input data is read, memory is allocated, data structures are initialized, and the communication scheme is defined. If the cost of this step is too high, any potential advantages gained in solving the problem may be diminished or even lost. Therefore, minimizing the cost of this stage is required to ensure overall efficiency. Before analyzing GLP ’s convergence, we compare the preprocessing time for GLP and LLP using the same cases discussed in Section 4.5. Figure 4.18 presents the ratio of preprocessing time for GLP relative to LLP. Unlike LLP, GLP requires additional operations, such as assembling the local linear system, computing contributions to the interface system, mapping these contributions, and handling communication. Therefore, preprocessing takes longer for GLP, with its time increasing as the number of domain partitions grows. The figure indicates that GLP ’s preprocessing time ranges between 1.34 and 4.33 times that of LLP . 106 4.7. Conclusions References [1] Hadrien Calmet, Alberto M. Gambaruto, Alister J. Bates, Mariano Vázquez, Guillaume Houzeaux and Denis J. Doorly. ‘Large-scale CFD simulations of the transitional and turbulent regime for the large human airways during rapid inhalation’. In: Computers in Biology and Medicine 69 (2016), pp. 166–180. issn: 0010-4825. doi:https://doi.org/10.1016/j. compbiomed.2015.12.003.url:https://www.sciencedirect.com/science/ article/pii/S0010482515003881. [2] R. de Olazábal, R. Borrell and O. Lehmkuhl. ‘An algebraic global linelet preconditioner for incompressible flow solvers’. In: Journal of Computational Physics 514 (2024), p. 113237. issn: 0021-9991. doi:https: //doi.org/10.1016/j.jcp.2024.113237.url:https://www.sciencedirect.com/ science/article/pii/S0021999124004868. [3] S. Gövert, D. Mira, M. Zavala-Ake, J.B.W. Kok, M. Vázquez and G. Houzeaux. ‘Heat loss prediction of a confined premixed jet flame using a conjugate heat transfer approach’. In: International Journal of Heat and Mass Transfer 107 (2017), pp. 882–894. issn: 0017-9310. doi:https: / / doi . org / 10 . 1016 / j . ijheatmasstransfer. 2016 . 10 . 122.url:https : //www.sciencedirect.com/science/article/pii/S0017931016317495. [4] O. Lehmkuhl, G. Houzeaux, H. Owen, G. Chrysokentis and I. Rodriguez. ‘A low-dissipation finite element scheme for scale resolving simulations of turbulent flows’. In: Journal of Computational Physics 390 (2019), pp. 51– 65. issn: 0021-9991. doi:https://doi.org/10.1016/j.jcp.2019.04.004.url: https://www.sciencedirect.com/science/article/pii/S0021999119302372. [5] Alfonso Santiago, Miguel Zavala-Aké, Ricard Borrell, Guillaume Houzeaux and Mariano Vázquez. ‘HPC compact quasi-Newton algorithm for interface problems’. In: Journal of Fluids and Structures 96 (2020), p. 103009. issn: 0889-9746. doi:https : / / doi . org / 10 . 1016 / j . jfluidstructs . 2020 . 103009.url:https : / / www. sciencedirect . com / science / article / pii / S088997461930934X. [6] Orlando Soto, Rainald Lohner and Fernando Camelli. ‘A Linelet preconditioner for incompressible flow solvers’. In: International Journal of Numerical Methods for Heat & Fluid Flow - INT J NUMER METHOD HEAT FL F 13 (Feb. 2003), pp. 133–147. doi:10.1108/ 09615530310456796. [7] Mariano Vázquez, Guillaume Houzeaux, Seid Koric, Antoni Artigues, Jazmin Aguado-Sierra, Ruth Arís, Daniel Mira, Hadrien Calmet, Fernando Cucchietti, Herbert Owen, Ahmed Taha, Evan Dering Burness, José María Cela and Mateo Valero. ‘Alya: Multiphysics engineering simulation toward exascale’. In: Journal of Computational Science 14 (2016). The Route to Exascale: Novel Mathematical Methods, Scalable Algorithms and Computational Science Skills, pp. 15–27. issn: 1877-7503. doi:https: //doi.org/10.1016/j.jocs.2015.12.007.url:https://www.sciencedirect. com/science/article/pii/S1877750315300521. 113 CHAPTER 5 Numerical Applications – Solving Real-World Scenarios Part of the contents of this chapter have been published as: R. de Olazábal, R. Borrell and O. Lehmkuhl. ‘An algebraic global linelet preconditioner for incompressible flow solvers’. In: Journal of Computational Physics 514 (2024), and disseminated in the following conferences: the 32nd Parallel Computational Fluid Dynamics Conference, Nice, France, 2021, the 15th JLESC Workshop, Bordeaux, France, 2023, the 1st Math 2 Product (M2P) conference for Emerging Technologies in Computational Science for Industry, Sustainability and Innovation, Taormina, Sicily, 2023, and the 23rd IACM Computational Fluids Conference, Santiago de Chile, 2025. In this chapter, we assess the performance of the GLP preconditioner in three large-scale computational fluid dynamics (CFD) simulations. These cases represent three real-world flow problems involving different geometries, Reynolds numbers, and turbulence modeling approaches. The objective is to evaluate how GLP performs in practical scenarios, particularly in terms of convergence behavior, computational cost, and parallel efficiency. The first test case investigates the flow over the 30P30N three-element high-lift airfoil, a well-established benchmark in aerodynamics. This simulation employs LES at a Reynolds number of Rec = 750 , 000 and involves a hybrid unstructured mesh with structured inflation layers around the airfoil surfaces. This configuration has been extensively analyzed in [6]. The second case is a DNS of flow through the Stanford diffuser at Re = 10 , 000, which serves as a canonical case for analyzing flow separation and turbulence in a confined domain. It was initially introduced by Cherry et al. [1] and further analyzed in [5]. 114 5.1. 30P30N Airfoil Finally, the third case examines the DrivAer model, a detailed representation of a realistic automobile, developed at the Institute of Aerodynamics and Fluid Mechanics at Technische Universität München (TUM) [2,3,4]. The simulation is performed at Re = 4 . 87 × 10 6 (based on vehicle length), using Wall-modelled Large Eddy Simulation (WMLES). The DrivAer case is the largest considered in this study, comprising 427 million elements and requiring up to 14,000 computational partitions. In all three cases, we employ an open integration rule. Given that GLP benefits significantly from a geometric approach (see, for example, Figure 1.4), we opted for a geometric-based linelet construction method in these simulations. The chapter is organized as follows. Sections 5.1, 5.2, and 5.3 present detailed analyses of the 30P30N airfoil, Stanford diffuser, and DrivAer model cases, respectively. Each section includes an evaluation of solver convergence, computational cost, and the impact of GLP on numerical stability. These analyses compare the performance of GLP against the LLP and the diagonal preconditioner. Following this, Section 5.4 provides a comprehensive performance analysis, comparing preprocessing and preconditioning costs across all three cases. Here, we highlight the communication overhead and load balance effects that influence overall computational efficiency. Finally, Section 5.5 presents a detailed scalability study, examining speed-up, communication efficiency, load balance, and parallel efficiency. This analysis extends the findings from Chapter 4, reinforcing key insights into how GLP scales across different problem sizes and partition counts. The results obtained in this chapter serve as a validation of the GLP method in large-scale CFD applications. While significant improvements in convergence are demonstrated, trade-offs related to preprocessing time and communication overhead are also identified. These findings provide guidance for future optimizations, particularly in hybrid preconditioning strategies and domain decomposition techniques. The chapter concludes with a discussion of these insights and their implications for the efficient solution of real-world fluid dynamics problems. 5.1 30P30N Airfoil Before analyzing overall convergence trends, we first examine the residual evolution during each iteration of the Preconditioned Conjugate Gradient (PCG) method for the high-lift wing case at the first time step. Figure 5.1 presents the residual convergence for the diagonal, LLP, and GLP preconditioners, with the latter two evaluated using both geometric and algebraic approaches. A key observation from Figure 5.1 is that while the algebraic GLP exhibits a significant improvement over the algebraic LLP, its advantage over the diagonal preconditioner diminishes after a residual of 3 × 10 −4 , beyond which the diagonal 115 5.1. 30P30N Airfoil Figure 5.1: Residual convergence over time for the PCG method applied to the 30P30N case at Rec = 750000. The figure compares the algebraic and geometric approaches for constructing the linelets. The simulation was performed on 1680 processes. preconditioner becomes the fastest. This behavior is consistent with Figures 4.9 and 4.10, which illustrate that the benefits of GLP over diagonal and LLP preconditioners become more pronounced when the proportion of elements within the boundary layer exceeds 40%. In this case, approximately 43% of the elements lie within the boundary layer, meaning the expected improvement is relatively moderate. This observation reinforces the findings from previous analyses—GLP is particularly effective in scenarios with a high percentage of elements in the boundary layer. Furthermore, this result demonstrates that the conclusions drawn from section remain valid even for complex geometries and higher Reynolds numbers. Another crucial point is the advantage of the geometric approach over the algebraic one. The geometric GLP converges significantly faster, particularly at lower tolerances. For instance, it reaches a residual tolerance of 10 −6 in 105.73 seconds, while the diagonal preconditioner, the second fastest method, requires 202.80 seconds. This translates into a speed-up of approximately 1.92 times, or, conversely, the geometric GLP completing in just 52% of the time required by the diagonal preconditioner. These findings align with Figures 1.4 and 1.5, further highlighting the superior efficiency of the geometric approach. 5.1.1 Computation and Communication Time Analysis To better understand the computational cost of the GLP preconditioner, we analyze the mean time per PCG iteration required for both computation and communication. This is compared against the LLP preconditioner as a function of the number of processes. Figure 5.2 presents the mean time per iteration for different process counts: 896, 1344, 1792, 2688, and 3584. These values 116 5.1. 30P30N Airfoil Figure 5.2: Mean time per iteration for the LLP and GLP preconditioners for the 30P30N case at Rec= 750000. represent the average over all threads. These results indicate that the communication step is the primary bottleneck in the GLP preconditioner. Meanwhile, the computation step takes at most twice the time the LLP preconditioner requires. This trend is further illustrated in Figure 5.4a, which shows the ratio of the time consumed by each step of the GLP compared to the corresponding mean time in the LLP approach. The increased cost of the GLP preconditioner aligns with the additional computational steps involved (see Section 2.2.2). Unlike LLP, the GLP approach explicitly solves the preconditioning system for interface nodes. Furthermore, in the GLP method, more one-node slices may be incorporated into the linelet system. These slices are simply treated with a diagonal preconditioner in the LLP approach. Figures 5.2 and 5.4a highlight that, on average, the computation time required by GLP is approximately twice that of LLP. Additionally, as expected, communication time increases with the number of processes, suggesting that parallel overhead becomes more pronounced at larger scales. 5.1.2 Computational Cost per PCG Iteration To further analyze the performance of the GLP preconditioner, we examine the time taken per iteration of the PCG method. Figure 5.3 presents the total time required for a single PCG iteration for both LLP and GLP. This data can be compared with Figure 5.4b, which illustrates the ratio of the total time taken by GLP relative to LLP. The results indicate that the GLP preconditioner is, in this case, at most 3.56 times slower than LLP. 117 5.1. 30P30N Airfoil Figure 5.3: Time per iteration for the LLP and GLP preconditioners for the 30P30N case at Rec= 750000. To gain further insight, Figure 5.4a provides a breakdown of the mean time per iteration, distinguishing between computation and communication times. It also includes the ratio of total time for GLP relative to LLP. These results are consistent with previous observations (see Figure 4.8). The plot in Figure 5.4a was computed by averaging the times over all processes, providing an overview of the general behavior. A similar pattern is observed in Figure 4.4, with the added distinction that computation and communication times are explicitly compared. However, this averaged perspective may not fully capture the total computational cost of GLP. The slowest process dictates the overall time in the presence of load imbalance. Given the nature of the linelet-based preconditioners (GLP and LLP), load imbalance can arise when linelets are concentrated in specific regions of the mesh, and the mesh partitioning does not account for the linelet distribution. This effect is illustrated in Figure 4.11. Future work can be done to consider load-balancing strategies in large-scale implementations. On the other hand, a similar analysis can be conducted for the preprocessing step. Figure 5.5 presents the average time required per process to complete the preprocessing stage for the LLP and GLP preconditioners. The results indicate that, on average, the computational cost of the preprocessing step in GLP is at most 1.95 times higher than that of LLP. However, the communication overhead significantly amplifies the overall cost, making the GLP preconditioner up to 7.35 times more expensive because of the communication time. However, when comparing the total time spent in the preprocessing stage, the load imbalance present in this stage makes the GLP at most 2.42 times more expensive than the required for the LLP. This suggests that while the additional communication steps required by GLP increase the preprocessing cost, their 118 5.1. 30P30N Airfoil (a) Quotient between the mean computation and communication times for LLP and GLP, alongside the total time ratio for GLP over LLP. (b) Total time quotient for the LLP and GLP preconditioners. Figure 5.4: Comparison of the mean time per iteration for the LLP and GLP preconditioners in the 30P30N case. 119 5.1. 30P30N Airfoil impact is partially mitigated by load-balancing effects across processes. After verifying that the behavior of the GLP preconditioner aligns with the results obtained in the previous sections, we now analyze its overall impact on the time required to solve the Navier-Stokes equations. Figure 5.6 presents the time per time step required to solve the Navier-Stokes equations for the 30P30N case at Rec = 750000 using the diagonal, LLP, and GLP preconditioners, considering 3584 partitions and tolerances of 10 2 and 10 6 . The results indicate that when employing the GLP, the convergence demonstrates the most stable performance, rapidly reaching a steady-state convergence time with minimal fluctuations. In contrast, the LLP preconditioner consistently requires more time per step and exhibits higher variability. On the other hand, the diagonal preconditioner shows large fluctuations in convergence time, particularly in the early iterations. This suggests that the GLP preconditioner accelerates convergence and enhances numerical stability in this specific scenario. Another way to visualize these results is to plot the cumulative time required to solve the Navier-Stokes equations. Figures 5.7a and 5.7b show the cumulative time required to perform n time steps for the cases presented in Figures 5.6a and 5.6b. These results confirm that the GLP preconditioner is the most efficient, followed by the LLP. A broader analysis can be performed by varying the number of processes and tolerances. Table 5.1 presents the ratio of the time required by GLP to that required by LLP to perform n time steps, where n is the last comparable time step. This concept is illustrated in Figure 5.8. Since the choice of preconditioner in the pressure correction equation influences the overall time required to solve the Navier-Stokes equations, the LLP and GLP cases may not complete the same number of time steps within the given simulation time frame. Thus, we take the last comparable time step reached by both preconditioners. For instance, in Figure 5.8, the LLP executes 209 time steps, while the GLP completes 303. Therefore, we compute the quotient between the time required by both cases to reach 209 time steps. For this scenario, the GLP takes only 29.6% of the time required by the LLP, demonstrating a substantial improvement. Table 5.1 further confirms that the most significant advantages of the GLP are obtained when using lower tolerances and a higher number of processes. This trend is consistent with the findings in Section 5.1 (see Figures 4.9 and 4.10). Blank spaces in the table indicate missing data. 120 5.1. 30P30N Airfoil (a) Quotient between the mean time spent in computation and communication for the preprocessing stage in the LLP and GLP preconditioners. (b) Total time quotient for the preprocessing stage in the GLP and LLP preconditioners. Figure 5.5: Comparison of the preprocessing stage for the LLP and GLP preconditioners in the 30P30N case. 121 5.1. 30P30N Airfoil (a) Comparison of the time per time step to solve the Navier-Stokes equations for the 30P30N case with a tolerance of 10 2 in the pressure correction equation and 3584 processes. (b) Same as Figure 5.6a but with a tolerance of 106. Figure 5.6: Comparison of the time per time step to solving the NavierStokes equations for the 30P30N case when using the diagonal, LLP, and GLP preconditioners. 122 5.3. DrivAer Model Figure 5.11: Time per time step required to solve the pressure correction equation for the DrivAer case using 14,000 processes and a tolerance of 10−5. Figure 5.12: Cumulative time required to solve the pressure correction equation for the DrivAer case with 14,000 processes and a tolerance of 10−5. 129 5.4. Performance Analysis Figure 5.13: Cumulative time required to solve the Navier-Stokes equations for the DrivAer case with 14,000 processes and a tolerance of 10−5. 5.4 Performance Analysis In this section, we evaluate the performance of the GLP preconditioner for the three studied cases, comparing it against the LLP approach, which serves as the benchmark. Since GLP is a generalization of LLP, any improvement observed in GLP should be measured relative to LLP. 5.4.1 Preprocessing Stage We begin by analyzing the preprocessing stage. Figure 5.14 presents the ratio between the preprocessing time required for GLP and LLP across all test cases. As expected, the preprocessing time scales with the problem size, with the DrivAer model—the largest case—exhibiting the longest preprocessing time. Notably, for this case, communication overhead significantly contributes to the total preprocessing time, leading to an increase of over an order of magnitude compared to LLP. To better understand this phenomenon, Figure 5.15 illustrates the mean preprocessing time per process for the DrivAer case. Here, the communication step is the dominant cost, confirming that communication overhead is the primary bottleneck in the preprocessing stage. On the other hand, Figure 5.16 presents the corresponding results for the Stanford diffuser case. Unlike the DrivAer model, the preprocessing time, in this case, is more balanced between computation and communication. This is expected, as the Stanford diffuser is the smallest case studied and involves fewer partitions. The smaller number of partitions results in less communication overhead, leading to a preprocessing time closer to that of LLP. These observations are consistent with the results discussed in Section 4.6 of Chapter 4, where we analyzed the preprocessing costs for a cavity flow case. In 130 5.4. Performance Analysis Figure 5.14: Comparison of preprocessing times for GLP and LLP across the three studied cases. Figure 5.15: Mean preprocessing time per process for the DrivAer case, showing the contribution of computation and communication. 131 5.4. Performance Analysis Figure 5.16: Mean preprocessing time per process for the Stanford diffuser case, highlighting the balanced cost between computation and communication. that analysis, we found that preprocessing time exhibited significant variability due to communication overhead, which was strongly influenced by the domain partitioning. The current results reinforce these findings: while computation time remains relatively stable across different cases, the communication cost is the main factor affecting preprocessing efficiency. As previously noted in Chapter 4, optimizing the domain partitioning strategy is one potential avenue for improving preprocessing performance. Since GLP is agnostic to mesh partitioning, an optimized partitioning approach that minimizes communication overhead while maintaining load balance could significantly enhance performance. This is especially relevant for large-scale cases like DrivAer, where communication costs dominate. Future work could explore partitioning strategies that take into account the spatial distribution of linelets, ensuring that linelets are not excessively fragmented across partitions. 5.4.2 Preconditioning Step In this section, we compare the computational cost of the preconditioning step for the three studied cases. Figure 5.17 presents the quotient of the total time spent in the preconditioning step for each case. It can be observed that the difference in computational cost between GLP and LLP is less pronounced in this step compared to the preprocessing stage. The overhead remains within the same order of magnitude as the LLP case, reinforcing the scalability of the GLP method. However, it is important to highlight that the most computationally expensive case is the 30P30N airfoil, which, paradoxically, exhibits the best convergence improvements. This is because a compromise is made between overcost in time and improvement in convergence. Depending on the case, one may outperform the other. This aligns with the results shown in Figure 5.1, demonstrating that the additional cost per iteration of the GLP does not necessarily correlate with inefficiency. 132 5.4. Performance Analysis Figure 5.17: Comparison of the mean time per process in the preconditioning step. Figure 5.18: Mean time per process in the computation and communication steps for the Stanford diffuser case. Figures 5.18 and 5.19 provide additional insights by breaking down the preconditioning step into computation and communication components for the Stanford diffuser and DrivAer cases, respectively. These figures show that in the diffuser and DrivAer cases, the communication step does not dominate the overall cost as much as in the 30P30N case. Despite this, in the former two cases, the GLP method does not significantly reduce iterations compared to LLP, suggesting that the underlying problem charac133 5.5. Scalability Analysis and Performance Metrics Figure 5.19: Mean time per process in the computation and communication steps for the DrivAer case. teristics play a crucial role in determining the effectiveness of the preconditioner. This behavior is consistent with the findings presented in Chapter 4, where the cavity flow case exhibited similar trends. In that case, the preconditioning step in GLP incurred an increased communication cost, mainly due to the communication between slices of linelets across different domain partitions. This further supports the idea that the primary source of overhead in GLP arises from interprocess communication, which is directly influenced by the domain decomposition strategy. Therefore, as in Chapter 4, the results reinforce the need for optimized partitioning strategies to better accommodate linelet structures and mitigate communication costs. Future work could explore adaptive partitioning techniques that account for linelet connectivity to further improve the efficiency of GLP in large-scale simulations. 5.5 Scalability Analysis and Performance Metrics To conclude this section, we evaluate the scalability and performance of the preconditioning step for the three studied cases: the DrivAer model, the 30P30N airfoil, and the Stanford diffuser. This analysis follows the methodology described in Chapter 4 (Section 4.5.2), where similar metrics were computed for a cavity flow case. 5.5.1 Speed Up Figure 5.20 and table 5.5 present the speed-up results for the three cases. The Stanford diffuser exhibits the closest behavior to ideal scaling, with both the 134 5.5. Scalability Analysis and Performance Metrics Figure 5.20: Speed-up for the DrivAer, 30P30N airfoil, and Stanford diffuser cases. GLP and LLP comparison. GLP and LLP preconditioners demonstrating strong scalability. Also, the GLP preconditioner has slightly worse scalability than the LLP, consistent with the fact that the GLP has a higher communication cost. In contrast, the DrivAer and 30P30N cases show more significant deviations from linear speed up, with the 30P30N airfoil performing the worst. This degradation aligns with previous observations (Figure 5.2), where the preconditioning step in the 30P30N case was found to be heavily dominated by communication time. Furthermore, the 30P30N airfoil exhibits the largest gap between the GLP and LLP cases, emphasizing the additional computational overhead introduced by the GLP preconditioner. For the DrivAer model, speed up remains reasonable but does not reach the levels observed in the diffuser case. This can be attributed to the sheer size of the problem (427 million elements) and the impact of load imbalance at higher process counts. 5.5.2 Communication Efficiency Figure 5.21 shows the communication efficiency for the three cases. It can be seen that the DrivAer model achieves the best communication efficiency, remaining above 0.9 even at high core counts. This indicates that within the tested number of partitions, the computational workload is well distributed and that communication overhead does not significantly hinder performance. In contrast, 135 5.5. Scalability Analysis and Performance Metrics Pr DrivAer 30P30N Diffuser GLP LLP GLP LLP GLP LLP 896 1.0 1.0 1.0 1.0 1344 1.0 1.0 1.20 1.20 1.42 1.47 1792 1.27 1.24 1.44 1.33 1.77 1.83 2688 1.77 1.79 1.81 1.93 2.53 2.64 3584 2.32 2.27 2.10 2.39 3.39 3.53 Table 5.5: Comparison of GLP and LLP Speed Up for different test cases and processor counts (Pr). Figure 5.21: Communication efficiency for the DrivAer, 30P30N airfoil, and Stanford diffuser cases. the 30P30N case experiences a noticeable decline in communication efficiency, dropping to 0.68 for the highest process counts. This trend is consistent with previous findings showing that the 30P30N airfoil spends the largest fraction of time in communication, which negatively affects its scalability. On the other hand, the Stanford diffuser case maintains a relatively stable communication efficiency of around 0.83 across all process counts, suggesting that its workload distribution is more balanced. The same results can be seen in detail in table 5.6. 5.5.3 Load Balance Figure 5.22 and table 5.7 present the load balance results for the three cases. The Stanford diffuser exhibits the best load balance, with the GLP and LLP configurations maintaining values above 0.55. In contrast, the 30P30N case 136 5.5. Scalability Analysis and Performance Metrics Pr DrivAer 30P30N Diffuser 896 0.81 0.85 1344 0.95 0.77 0.85 1792 0.95 0.85 0.83 2688 0.95 0.73 0.83 3584 0.97 0.68 0.85 Table 5.6: Communication Efficiency comparison of DrivAer, 30P30N, and Diffuser cases across different processor counts (Pr). experiences the worst load imbalance, with values dropping to 0.32 as the number of processes increases. This is likely due to the highly anisotropic nature of its mesh, which makes it more challenging to distribute computational workload evenly across processes. Interestingly, despite achieving the highest communication efficiency, the DrivAer model also suffers from a noticeable load imbalance, with values fluctuating between 0.40 and 0.51. This suggests that certain processes are underutilized due to workload variations. The observed imbalance is consistent with the behavior shown in Figure 4.11, where linelets concentrated in specific regions of the domain lead to increased computational costs for the processes handling those sections. However, the GLP preconditioner helps mitigate this issue by redistributing some of the computational workload associated with interface nodes across all processes, leading to improved balance. This effect is reflected in Figure 5.22, where the GLP configuration achieves better load balance than LLP for both the Stanford diffuser and DrivAer cases (see also Figure 4.14b). The 30P30N case, however, shows only a marginal improvement in load balance with GLP compared to LLP, indicating that the preconditioner’s balancing effect is less pronounced in this scenario. Pr DrivAer 30P30N Diffuser GLP LLP GLP LLP GLP LLP 896 0.50 0.55 0.65 0.60 1344 0.51 0.46 0.42 0.44 0.63 0.60 1792 0.48 0.44 0.34 0.37 0.59 0.56 2688 0.45 0.42 0.34 0.35 0.57 0.54 3584 0.44 0.40 0.32 0.33 0.55 0.54 Table 5.7: Load balance comparison of GLP and LLP across different processor counts (Pr) for DrivAer, 30P30N, and Diffuser test cases. 5.5.4 Parallel Efficiency Parallel efficiency, shown in Figure 5.23 and in Table 5.8, reflects the combined effects of speed-up, communication efficiency, and load balance, providing a comprehensive measure of scalability. The Stanford diffuser achieves the highest parallel efficiency among the three cases, with both GLP and LLP 137 5.5. Scalability Analysis and Performance Metrics Figure 5.22: Load balance for the DrivAer, 30P30N airfoil, and Stanford diffuser cases. GLP and LLP comparison. configurations maintaining values above 0.47. The DrivAer model follows, though its efficiency declines as the number of processes increases. As expected, the 30P30N case exhibits the lowest parallel efficiency, dropping below 0.22 at the highest core counts. These results align with the trends observed in Figure 5.2, highlighting the substantial communication overhead incurred by the 30P30N case. In all three cases, the GLP configuration consistently demonstrates lower parallel efficiency than the LLP. This confirms previous findings that while the GLP significantly accelerates convergence, its added computational and communication overhead reduces overall efficiency. Pr DrivAer 30P30N Diffuser GLP LLP GLP LLP GLP LLP 896 0.40 0.55 0.55 0.60 1344 0.48 0.46 0.33 0.44 0.53 0.60 1792 0.46 0.44 0.29 0.37 0.49 0.56 2688 0.43 0.42 0.25 0.35 0.47 0.54 3584 0.42 0.40 0.22 0.33 0.47 0.54 Table 5.8: Parallel efficiency comparison of GLP and LLP across different processor counts (Pr) for the DrivAer, 30P30N, and Diffuser test cases. 138