Parallel Solver Research Articles

An open-source, adaptive solver for particle-resolved simulations with both subcycling and non-subcycling methods

We present the IAMReX (incompressible flow with adaptive mesh refinement for the eXascale), an adaptive and parallel solver for particle-resolved simulations on the multi-level grid. The fluid equations are solved using a finite-volume scheme on the block-structured semi-staggered grids with both subcycling and non-subcycling methods. The particle-fluid interaction is resolved using the multidirect forcing immersed boundary method. The associated Lagrangian markers used to resolve fluid-particle interface only exist on the finest-level grid, which greatly reduces memory usage. The volume integrals are numerically calculated to capture the free motion of particles accurately, and the repulsive potential model is also included to account for the particle–particle collision. We demonstrate the versatility, accuracy, and efficiency of the present multi-level framework by simulating fluid-particle interaction problems with various types of kinematic constraints. The cluster of monodisperse particles case is presented at the end to show the capability of the current solver in handling multiple particles. It is demonstrated that the three-level AMR (Adaptive Mesh Refinement) simulation leads to a 72.46% grid reduction compared with the single-level simulation. The source code and testing cases used in this work can be accessed at https://github.com/ruohai0925/IAMR/tree/development. Input scripts and raw postprocessing data are also available for reproducing all results.

NUMA-aware parallel sparse LU factorization for SPICE-based circuit simulators on ARM multi-core processors

In circuit simulators that resemble the Simulation Program with Integrated Circuit Emphasis (SPICE), one of the most crucial steps is the solution of numerous sparse linear equations generated by frequency domain analysis or time domain analysis. The sparse direct solvers based on lower-upper (LU) factorization are extremely time-consuming, so their performance has become a significant bottleneck. Despite the existence of some parallel sparse direct solvers for circuit simulation problems, they remain challenging to adapt in terms of performance and scalability in the face of rapidly evolving parallel computers with multiple NUMA hardware based on ARM architecture. In this paper, we introduce a parallel sparse direct solver named HLU, which re-examines the performance of the parallel algorithm from the viewpoint of parallelism in pipeline mode and the computing efficiency of each task. To maximize task-level parallelism and further minimize the thread waiting time, HLU devises a fine-grained scheduling method based on an elimination tree in pipeline mode, which employs depth-first search (DFS-like) to iteratively search for parent tasks and then place dependent tasks in the same task queue. HLU also suggests two NUMA node affinity strategies: thread affinity optimization based on NUMA nodes topology to guarantee computational load balancing and data affinity optimization to enable effective memory placement when threads access data. The rationality and effectiveness of the sparse solver HLU are validated by the SuiteSparse Matrix Collection. In comparison with KLU and NICSLU, the experimental results and analysis show that HLU attains a speedup of up to 9.14× and 1.26x (geometric mean) on a Huawei Kunpeng 920 Server, respectively.

Comprehensive Evaluation of the Massively Parallel Direct Simulation Monte Carlo Kernel “Stochastic Parallel Rarefied-Gas Time-Accurate Analyzer” in Rarefied Hypersonic Flows—Part A: Fundamentals

The Direct Simulation Monte Carlo (DSMC) method, introduced by Graeme Bird over five decades ago, has become a crucial statistical particle-based technique for simulating low-density gas flows. Its widespread acceptance stems from rigorous validation against experimental data. This study focuses on four validation test cases known for their complex shock–boundary and shock–shock interactions: (a) a flat plate in hypersonic flow, (b) a Mach 20.2 flow over a 70-degree interplanetary probe, (c) a hypersonic flow around a flared cylinder, and (d) a hypersonic flow around a biconic. Part A of this paper covers the first two cases, while Part B will discuss the remaining cases. These scenarios have been extensively used by researchers to validate prominent parallel DSMC solvers, due to the challenging nature of the flow features involved. The validation requires meticulous selection of simulation parameters, including particle count, grid density, and time steps. This work evaluates the SPARTA (Stochastic Parallel Rarefied-gas Time-Accurate Analyzer) kernel’s accuracy against these test cases, highlighting its parallel processing capability via domain decomposition and MPI communication. This method promises substantial improvements in computational efficiency and accuracy for complex hypersonic vehicle simulations.

A parallel solver framework for fully implicit monolithic fluid-structure interaction

A parallel solver framework for fully implicit monolithic fluid-structure interaction

Exploring the impact of adenoid obstruction on the human nasal airflow through parallel computational fluid dynamics

AbstractAdenoid obstruction is a prevalent condition that can significantly alter nasal airflow patterns, leading to various respiratory symptoms. This study aims to explore the influence of adenoid obstruction on human nasal airflow using computational fluid dynamics (CFD) simulations. We developed a realistic nasal cavity model to simulate both steady‐state inhalation and transient breathing cycles, accounting for various levels of obstruction. The three‐dimensional, incompressible Navier‐Stokes equations were discretized employing a stabilized finite element method on unstructured meshes, along with an implicit second‐order backward differentiation formula. The numerical solution was obtained using a Newton‐Krylov‐Schwarz algorithm‐based parallel solver. Our findings indicate that obstruction of the adenoid cross‐sectional area below 50% minimally affects nasal airflow, whereas noticeable changes occurred as obstructions exceeded the 50% threshold. Complete blockage was observed at levels greater than 70%, and inhalation difficulties could arise even at the 50% obstruction mark during stimulation or extreme inhalation scenarios. These results substantiate the clinical guidelines that advocate a 50% obstruction ratio as the demarcation between conservative treatment and surgery. Furthermore, this study offers valuable insights into the intricate relationship between the adenoid obstruction and nasal airflow, potentially enhancing clinical decision‐making in the diagnosis and treatment of nasal obstruction.

Efficient parallel solver for rarefied gas flow using GSIS

Recently, the general synthetic iterative scheme (GSIS) has been proposed to find the steady-state solution of the Boltzmann equation in the whole range of gas rarefaction, where its fast-converging and asymptotic-preserving properties lead to the significant reduction of iteration numbers and spatial cells in the near-continuum flow regime. However, the efficiency and accuracy of GSIS have only been demonstrated in two-dimensional problems with small numbers of spatial cells and discrete velocities. Here, a large-scale parallel computing strategy is designed to extend the GSIS to three-dimensional flow problems, including the supersonic flows which are usually difficult to solve by the discrete velocity method. Since the GSIS involves the calculation of the mesoscopic kinetic equation which is defined in six-dimensional phase-space, and the macroscopic high-temperature Navier–Stokes–Fourier equations in three-dimensional physical space, the proper partition of the spatial and velocity spaces, and the allocation of CPU cores to the mesoscopic and macroscopic solvers, are the keys to improving the overall computational efficiency. These factors are systematically tested to achieve optimal performance, up to 100 billion spatial and velocity grids. For hypersonic flows around the Apollo reentry capsule, the X38-like vehicle, and the space station, our parallel solver can obtain the converged solution within one hour.

Parameter estimation in a whole-brain network model of epilepsy: Comparison of parallel global optimization solvers.

The Virtual Epileptic Patient (VEP) refers to a computer-based representation of a patient with epilepsy that combines personalized anatomical data with dynamical models of abnormal brain activities. It is capable of generating spatio-temporal seizure patterns that resemble those recorded with invasive methods such as stereoelectro EEG data, allowing for the evaluation of clinical hypotheses before planning surgery. This study highlights the effectiveness of calibrating VEP models using a global optimization approach. The approach utilizes SaCeSS, a cooperative metaheuristic algorithm capable of parallel computation, to yield high-quality solutions without requiring excessive computational time. Through extensive benchmarking on synthetic data, our proposal successfully solved a set of different configurations of VEP models, demonstrating better scalability and superior performance against other parallel solvers. These results were further enhanced using a Bayesian optimization framework for hyperparameter tuning, with significant gains in terms of both accuracy and computational cost. Additionally, we added a scalable uncertainty quantification phase after model calibration, and used it to assess the variability in estimated parameters across different problems. Overall, this study has the potential to improve the estimation of pathological brain areas in drug-resistant epilepsy, thereby to inform the clinical decision-making process.

A coupled fully implicit method for solving drag-adjoint equations of steady incompressible flows

Numerical instability is one of the main concerns when applying adjoint-based gradient computation in aerodynamic design optimization for a variety of engineering problems. An adjoint pressure-adjoint velocity-coupled fully implicit method for solving the adjoint equations of steady incompressible flows is developed in the present work, taking the aerodynamic drag as the optimization target. The finite volume method for cell-centered, collocated variables on unstructured grids is used for spatial discretization. A stabilization term based on pressure-weighted interpolation (PWI) of the adjoint velocity is added to the adjoint continuity equation to avoid spurious oscillations of adjoint pressure. All the terms of the continuous adjoint equations are discretized implicitly, including the adjoint transpose convection (ATC) term, which is regarded as one of the major sources of numerical instability. The resulting discretized equations are integrated into one linear system, and the adjoint fields can be obtained simultaneously. To prevent the diagonal element of the coefficient matrix of the linear system from getting too small, a singular value decomposition is applied to the diagonal block of the coefficient matrix, and a limiter is applied to the diagonal elements. The condition number of the coefficient matrix may be quite large, so a parallel sparse linear solver that effectively combines direct and iterative methods through a Schur complement approach is applied to solve the linear system. Compared to the staggered iterative solver, no damping or masking is needed for the ATC term in the proposed fully implicit method. There is no need to sacrifice numerical accuracy for numerical stability. This method has been tested in several benchmark cases, and the numerical results show good numerical stability and accuracy. The present method provides a sound basis for industrial applications of adjoint-based aerodynamic optimization methods.

GPU-accelerated approach for 2D fracture analysis of structures combining finite particle method and cohesive zone model

The fracture analysis of structures typically involves strong discontinuity and nonlinear behaviors, and it requires fine meshing and small time steps, making it time-consuming. This study proposes an approach accelerated by graphics processing units (GPU) for two-dimensional (2D) fracture analysis of structures combining finite particle method (FPM) and cohesive zone model (CZM). Specifically, the equations of motion of particles are presented, and the internal forces of planar triangular elements are derived based on FPM. Subsequently, the internal forces of four-particle cohesive element are derived to describe the CZM, and the linear elastic stage, softening stage, and failure stage are depicted according to the traction-separation criteria. An explicit GPU-based FPM analysis strategy for structural fracture analysis is then developed, and the parallel solvers for particles, triangular elements, and cohesive elements are implemented. A quasi-static fracture analysis of a concrete beam and a dynamic fracture analysis of the Kalthoff problem is conducted to verify the effectiveness of the proposed approach. The obtained crack propagation path, load vs. displacement curves, and crack propagation speed are in good agreement with the experimental results and the results in the literature. The efficiency of the proposed approach is also investigated. The achieved maximum speedup ratio of the GPU-accelerated approach to the central processing unit (CPU)-based approach reaches around 24, and the achieved speedup ratio relative to the commercial finite element software Abaqus/Explicit reaches around 11.5, demonstrating the computational efficiency of the proposed approach.

An integrated three-dimensional aeromechanical analysis for the prediction of stresses on modern coaxial rotors

This paper presents the first application of an Integrated Three-Dimensional aeromechanical analysis—defined as the coupled solution of three-dimensional finite element-based structural dynamics with a three-dimensional Reynolds-Averaged Navier-Stokes-based fluid dynamics—to predict the stresses on a modern coaxial rotor. The coupling was carried out with the University of Maryland’s structural dynamic solver X3D and the U.S. Army’s CREATETM–AV Helios suite of fluid dynamic solvers. A modern four-bladed hingeless coaxial rotor model—inspired by the gross dimensions of the Sikorsky S-97 Raider but generic and open source otherwise—is developed as a demonstration test case. The new structural solver is driven by parallel and scalable solvers and advanced high performance computing. It is enabled by high-order three-dimensional brick finite elements unified with multibody dynamics, integrated aerodynamics, and a special 3D-to-1D fluid-structure interfaces refines the power of delta-coupling procedure while retaining the advantages of existing CFD mesh motion schemes. The analysis predicts the three-dimensional stresses on the rotor blades and hub, together with the deformations, airloads, and wake, in an integrated manner. Two flight conditions are studied: a low-speed flight at 37 knots and a high-speed flight at 150 knots. Interesting three-dimensional unsteady stress patterns are revealed all across the blade but particularly inboard of 50% rotor radius—patterns that change from flight to flight and have remained invisible until now—since they could neither be predicted nor measured in flight. The maximum axial stresses exhibited 3/rev variation at low speed, and 2/rev variation at high speed flight. The lower rotor carried higher oscillatory stress burden at low speed, whereas both the rotors shared the same stress burden at high speed flight. The key conclusion is that such analysis is now indeed possible, and the stress patterns they reveal provide deeper insights into the dynamics of advanced rotors, and these might provide a path toward mitigating them in the future.

Graph and convolutional neural network coupling with a high-performance large-eddy simulation solver

Computational Fluid Dynamics (CFD) traditionally relies on long-standing numerical simulation strategies for the Navier–Stokes equations. Recently, interest in data-driven hybrid CFD solvers has spiked, leveraging pre-computed datasets to enhance various weak links inside existing solvers, such as closure models, under-resolved physics, or even to guide numerical resolution strategies. Running these hybrid solvers, notably in High Performance Computing (HPC) environments, presents specific challenges. In particular, context-aware deep learning (e.g. Convolutional (CNN) and Graph (GNN) Neural Networks) is promising for this task, but requires passing data representations between the physics solver and the neural network. In relevant industrial configurations, CFD meshes can be Cartesian but highly irregular, or unstructured, both of which do not match the pixel/voxel structure needed to run CNNs. In addition, discrepancies in programming language and libraries are common between CFD and machine learning applications. This work explores the many challenges of running a parallel hybrid solver in an HPC context, through the coupling of the AVBP CFD solver with neural networks in turbulent combustion and wall friction modeling applications. The knowledge gained is showcased in this article, as well as assembled in an actionable open-source library.

On cost-efficient parallel iterative solvers for 3D frequency-domain seismic multisource viscoelastic anisotropic wave modeling

Solving large sparse linear systems in 3D frequency-domain seismic wave modeling, especially in viscoelastic anisotropic media, poses significant challenges due to the increasing number of discrete moduli and nonzero elements in the linear system matrix. The computational load surpasses that of acoustic or viscoacoustic media, making it even more challenging when dealing with multisource problems. Popular scientific tools for solving a linear system, such as multifrontal massively parallel sparse direct solver (MUMPS), structured matrix package (STRUMPACK), and portable extensible toolkit for scientific computation (PETSc), can be used, but their applicability to our specific problem has not been comprehensively evaluated. Our study aims to tackle the challenges in solving large, sparse, complex-valued symmetric linear systems with multiple right-hand side vectors for 3D frequency-domain seismic wave modeling. We have adopted the preconditioned conjugate gradient iterative algorithms as the foundation for our research, introducing two highly cost-effective parallel iterative solvers: the parallel symmetric successive overrelaxation conjugate gradient (P-SSORCG) and the parallel incomplete Cholesky conjugate gradient (P-ICCG). These novel solvers are subjected to a comprehensive comparative analysis against well-established scientific tools, such as MUMPS, STRUMPACK, and PETSc, in the context of 3D frequency-domain seismic wave modeling. We show their promising performances in a practical 3D SEG/EAGE overthrust model and demonstrate that the grouped P-SSORCG offers an efficient alternative to parallel direct solvers, particularly in situations wherein computational resources are limited.

Integrated Three-Dimensional Airloads and Stresses on Lift-Offset Coaxial Rotors at Extreme Speeds

This paper investigates the three-dimensional (3-D) dynamic stresses on a modern four-bladed hingeless coaxial rotor—inspired by the gross dimensions of the Sikorsky S-97 Raider at extreme flight speeds. The stresses are obtained using integrated 3-D (I3D) aeromechanical analysis—defined as the coupling of 3-D finite element-based structural dynamics with 3-D Reynolds-averaged Navier–Stokes–based fluid dynamics. The coupling was carried out with the University of Maryland/U.S. Army structural dynamic solver X3D and the U.S. Army CREATETM–AV Helios suite of fluid dynamic solvers. The new structural analysis is both enabled and driven by advanced high-performance computing, parallel and scalable solvers, high-order 3-D brick finite elements unified with multibody dynamics, integrated aeromechanics, and a special 3D-to-1D fluid–structure interface that refines the power of the delta-coupling procedure while retaining the advantages of existing computational fluid dynamics mesh motion schemes. The analysis is carried out at 220 knots (μ=0.5)—the cruise speed of the S-97 Raider without reduced tip speed—in order to study the stresses in extreme conditions. At such high speeds, the blade lift is dominated by the complex tip vortex roll-up, and the pitching moments and drag are dominated by the unsteady transonic shocks at the tip. Interesting 3-D dynamic stress patterns are revealed all across the blade that have remained invisible until now since they could neither be predicted nor measured in flight. The key conclusion is that such high-fidelity analysis is now indeed possible and, in fact, necessary to get deeper insights into the dynamics of coaxial rotors at extreme speeds.

Parallel interior-point solver for block-structured nonlinear programs on SIMD/GPU architectures

ABSTRACT We investigate how to port the standard interior-point method to new exascale architectures for block-structured nonlinear programs with state equations. Computationally, we decompose the interior-point algorithm into two successive operations: the evaluation of the derivatives and the solution of the associated Karush-Kuhn-Tucker (KKT) linear system. Our method accelerates both operations using two levels of parallelism. First, we distribute the computations on multiple processes using coarse parallelism. Second, each process uses SIMD/GPU accelerators locally to accelerate the operations using fine-grained parallelism. The KKT system is reduced by eliminating the inequalities and the state variables from the corresponding equations. We demonstrate our method's capability on the supercomputer Polaris, a testbed for the future exascale Aurora system. Each node is equipped with four GPUs, a setup amenable to our two-level approach. Our experiments on the stochastic optimal power flow problem show that the reduction method is 50x faster than the sparse linear solver HSL MA57 running in serial on the CPU, and 6x faster than Pardiso running in parallel on CPU on the same number of processes.

Development and validation of a new one-dimensional parallel transient system solver with coordinate transformation method developed in OpenFOAM

Passive containment cooling system is an important component of both 3rd and 4th generation nuclear power plants. To simulate two-phase flow in loop-type passive containment cooling system, a new one-dimensional parallel transient system solver, ComspaFOAM-SYS1D, was developed based on OpenFOAM. Compared to the previously developed one-dimensional solver, the new solver can be accelerated using multithreading. Additionally, a new coordinate transformation method was used in ComspaFOAM-SYS1D, enabling OpenFOAM to simulate a one-dimensional loop with bends. Faucet problem was used to verify the basic equations and models in the new solver. Validation was undertaken based on the GENEVA facility and the passive containment cooling system comprehensive performance test facility. In the assessments, representative parameters such as void fraction, pressure, temperature, phase velocity, and oscillation period were compared. Regarding efficiency, the simulation results indicated that the parallel speedup of the new solver could approach a value of approximately 2.0. In terms of accuracy, the relative deviations of time-averaged and transient extremum values of all the simulation results from the experimental results were within 15.0%.

An offline coupling approach for efficient SPH simulations of long-distance tsunami events using wave source boundary condition

The entire procedure of tsunami events consists of wave generation and propagation toward the coastal structures encompassing the long-range system domain. Therefore, instead of an integrated continuous simulation, two or several consecutive simulations are preferred for numerical efficiency, where the hydrodynamic information within a particular domain from the preceding simulation can be utilized for the wave sources in the subsequent simulation. In this study, as an offline coupling approach, the wave source boundary condition is proposed based on the particle-based Lagrangian framework. This boundary condition is integrated with the solution algorithm of the smoothed particle hydrodynamics (SPH) in-house parallel code. In the dambreak simulation, it is validated that the wave source boundary condition can encompass both the forward and reverse flows, where the numerical results of the full and partial simulations show good agreement in terms of the wave profile, velocity distributions, and measured pressure histories. Thereafter, the full and partial simulations of a tsunami experiment are carried out. The numerical results in the partial simulation agree well with those obtained in the full simulation in terms of the particle distributions, physical quantities, and loading histories (e.g., pressure, acceleration) on the shoreline structure. It is demonstrated that the proposed offline coupling approach using the wave source boundary condition is compatibly utilized with the developed in-house SPH parallel solver, facilitating the detailed FSI simulations within a reduced domain of interest for long-distance tsunami-like phenomena.

Toward HPC application portability via C++ PSTL: the Gaia AVU-GSR code assessment

The computing capacity needed to process the data generated in modern scientific experiments is approaching ExaFLOPs. Currently, achieving such performances is only feasible through GPU-accelerated supercomputers. Different languages were developed to program GPUs at different levels of abstraction. Typically, the more abstract the languages, the more portable they are across different GPUs. However, the less abstract and co-designed with the hardware, the more room for code optimization and, eventually, the more performance. In the HPC context, portability and performance are a fairly traditional dichotomy. The current C++ Parallel Standard Template Library (PSTL) has the potential to go beyond this dichotomy. In this work, we analyze the main performance benefits and limitations of PSTL using as a use-case the Gaia Astrometric Verification Unit-Global Sphere Reconstruction parallel solver developed by the European Space Agency Gaia mission. The code aims to find the astrometric parameters of ∼108\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\sim10^8$$\\end{document} stars in the Milky Way by iteratively solving a linear system of equations with the LSQR algorithm, originally GPU-ported with the CUDA language. We show that the performance obtained with the PSTL version, which is intrinsically more portable than CUDA, is comparable to the CUDA one on NVIDIA GPU architecture.

Uncertainty quantification in the Henry problem using the multilevel Monte Carlo method

We investigate the applicability of the well-known multilevel Monte Carlo (MLMC) method to the class of density-driven flow problems, in particular the problem of salinisation of coastal aquifers. As a test case, we solve the uncertain Henry saltwater intrusion problem. Unknown porosity, permeability and recharge parameters are modelled by using random fields. The classical deterministic Henry problem is non-linear and time-dependent, and can easily take several hours of computing time. Uncertain settings require the solution of multiple realisations of the deterministic problem, and the total computational cost increases drastically. Instead of computing of hundreds random realisations, typically the mean value and the variance are computed. The standard methods such as the Monte Carlo or surrogate-based methods are a good choice, but they compute all stochastic realisations on the same, often, very fine mesh. They also do not balance the stochastic and discretisation errors. These facts motivated us to apply the MLMC method. We demonstrate that by solving the Henry problem on multi-level spatial and temporal meshes, the MLMC method reduces the overall computational and storage costs. To reduce the computing cost further, parallelization is performed in both physical and stochastic spaces. To solve each deterministic scenario, we run the parallel multigrid solver ug4 in a black-box fashion.

A parallel solver for fluid–structure interaction problems with Lagrange multiplier

The aim of this work is to present a parallel solver for a formulation of fluid–structure interaction (FSI) problems which makes use of a distributed Lagrange multiplier in the spirit of the fictitious domain method. The fluid subproblem, consisting of the non-stationary Stokes equations, is discretized in space by Q2–P1 finite elements, whereas the structure subproblem, consisting of the linear or finite incompressible elasticity equations, is discretized in space by Q1 finite elements. A first order semi-implicit finite difference scheme is employed for time discretization. The resulting linear system at each time step is solved by a parallel GMRES solver, accelerated by block diagonal or triangular preconditioners. The parallel implementation is based on the PETSc library. Several numerical tests have been performed on Linux clusters to investigate the effectiveness of the proposed FSI solver.

Massively Parallel Hybrid TLM-PEEC Solver and Model Order Reduction for 3D Nonlinear Electromagnetic Transient Analysis

Massively Parallel Hybrid TLM-PEEC Solver and Model Order Reduction for 3D Nonlinear Electromagnetic Transient Analysis