Scalable Solver Research Articles

Kernel multigrid on manifolds

Kernel methods for solving partial differential equations work coordinate-free on the surface and yield high approximation rates for smooth solutions. Localized Lagrange bases have proven to alleviate the computational complexity of usual kernel methods for data fitting problems, but the efficient numerical solution of the ill-conditioned linear systems of equations arising from kernel-based Galerkin solutions to PDEs is a challenging problem which has not been addressed in the literature so far. In this article we apply the framework of the geometric multigrid method with a τ≥2-cycle to scattered, quasi-uniform point clouds on the surface. We show that the resulting solver can be accelerated by using the Lagrange function decay and obtain satisfying convergence rates by a rigorous analysis. In particular, we show that the computational cost of the linear solver scales log-linear in the degrees of freedom.

Finite elements for Matérn-type random fields: Uncertainty in computational mechanics and design optimization

This work highlights an approach for incorporating realistic uncertainties into scientific computing workflows based on finite elements, focusing on prevalent applications in computational mechanics and design optimization. We leverage Matérn-type Gaussian random fields (GRFs) generated using the SPDE method to model aleatoric uncertainties, including environmental influences, variating material properties, and geometric ambiguities. Our focus lies on delivering practical GRF realizations that accurately capture imperfections and variations and understanding how they impact the predictions of computational models as well as the shape and topology of optimized designs. We describe a numerical algorithm based on solving a generalized SPDE to sample GRFs on arbitrary meshed domains. The algorithm leverages established techniques and integrates seamlessly with the open-source finite element library MFEM and associated scientific computing workflows, like those found in industrial and national laboratory settings. Our solver scales efficiently for large-scale problems and supports various domain types, including surfaces and embedded manifolds. We showcase its versatility through biomechanics and topology optimization applications, emphasizing the potential to influence these domains. The flexibility and efficiency of SPDE-based GRF generation empowers us to run large-scale optimization problems on 2D and 3D domains, including finding optimized designs on embedded surfaces, and to generate design features and topologies beyond the reach of conventional techniques. Moreover, these capabilities allow us to model and quantify geometric uncertainties on reconstructed submanifolds, such as the interpolated surfaces of cerebral aneurysms provided by postprocessing CT scans. In addition to offering benefits in these specific domains, the proposed techniques transcend specific applications and generalize to arbitrary forward and backward problems in uncertainty quantification involving finite elements.

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.

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.

Efficient multimode vectorial nonlinear propagation solver beyond the weak guidance approximation

In this paper, we present an efficient numerical model able to solve the vectorial nonlinear pulse propagation equation in circularly symmetric multimode waveguides. The algorithm takes advantage of the conservation of total angular momentum of light upon propagation and takes into account the vectorial nature of the propagating modes, making it particularly relevant for studies in ring-core fibers. While conventional propagation solvers exhibit a computational complexity scaling as Nmode4, where Nmode is the number of considered modes, the present solver scales as Nmode3/2. As a first example, it is shown that orbital angular momentum modulation instability processes take place in ring-core fibers in realistic conditions. Finally, it is predicted that the modulation instability process is followed by the appearance of breather-like angular structures.

Salome_cfd: An HPC framework for thermalhydraulics and multiphysics powered by the open-source Computational Fluid Dynamics Toolbox code_saturne

EDF R&D lab Chatou, Fluid Mechanics Energy and Environment department (MFEE), has been developping an HPC oriented Computational Fluid Dynamics toolbox since the late 1990s. This toolbox is based on the open-source Computational Fluid Dynamics (CFD) and multiphysics solver code_saturne. It has a proprietary extension, neptune_cfd, which is a Computational Multiphase Fluid Dynamics solver dealing with multi-phase flows. After a brief overview of code_saturne and neptune_cfd architecture and HPC capabilites, the present paper will focus on the effort made at EDF R&D to harness HPC clusters for thermalhydraulics at an industrial scale. This is made available by codes built from the start with massively-parallel architectures in mind, thus leading to highly scalable solvers which are able to handle tens of billions of cells. User-friendly access and interoperability are pivotal for engineers to take advantage of these HPC capabilities. Among common applications requiring HPC resources, a focus is first made on frontier scale simulations, then on parametric studies and uncertainty quantifications experiments, where a twolevel parallelism is needed : a parallel Design of Experiment evaluation drives computations themselves (massively) parallel.

Low Mach number aeroacoustic computations of airfoils in application with Pierce's wave equation

The airflows generated by drones are fairly different from those arising for conventional aeronautical architectures. They are low Mach number, relatively low Reynolds number and transient regimes play a dominant role. The goals of this work are twofold. This study firstly seeks to present an aeroacoustic model built on Pierce's wave equation, which is adapted to the specificity of this physic, and which is consistent and robust. The second part of this contribution then focuses on the use and the validation of a high-order highly scalable solver relying on an open source C++ library to solve the large scales of the flow turbulence by means of an explicit filtering. The flow corresponding to a benchmarked Re = 5 · 104 airfoil in full stall is computed and results are briefly discussed.

Coupling Conduction, Convection and Radiative Transfer in a Single Path-Space: Application to Infrared Rendering

In the past decades, Monte Carlo methods have shown their ability to solve PDEs, independently of the dimensionality of the integration domain and for different use-cases (e.g. light transport, geometry processing, physics simulation). Specifically, the path-space formulation of transport equations is a key ingredient to define tractable and scalable solvers, and we observe nowadays a strong interest in the definition of simulation systems based on Monte Carlo algorithms. We also observe that, when simulating combined physics (e.g. thermal rendering from a heat transfer simulation), there is a lack of coupled Monte Carlo algorithms allowing to solve all the physics at once, in the same path space, rather than combining several independent MC estimators, a combination that would make the global solver critically sensitive to the complexity of each simulation space. This brings to our proposal: a coupled, single path-space, Monte Carlo algorithm for efficient multi-physics problems solving. In this work, we combine our understanding and knowledge of Physics and Computer Graphics to demonstrate how to formulate and arrange different simulation spaces into a single path space. We define a tractable formalism for coupled heat transfer simulation using Monte Carlo, and we leverage the path-space construction to interactively compute multiple simulations with different conditions in the same scene, in terms of boundary conditions and observation time. We validate our proposal in the context of infrared rendering with different thermal simulation scenarios: e.g., room temperature simulation, visualization of heat paths within materials (detection of thermal bridges), heat diffusion capacity of thermal exchanger. We expect that our theoretical framework will foster collaboration and multidisciplinary studies. The perspectives this framework opens are detailed and we suggest a research agenda towards the resolution of coupled PDEs at the interface of Physics and Computer Graphics.

PaScaL_TCS: A versatile solver for large-scale turbulent convective heat transfer problems with temperature-dependent fluid properties

PaScaL_TCS is an efficient and scalable solver for large-scale direct numerical simulations of natural convective flow that considers temperature-dependent fluid properties. For increased scalability on a massive-scale distributed memory system, the solver decomposes the computational domain into a cubic sub-domain. The numerical procedure is based on the monolithic projection method with staggered time discretization [1], which is an implicit and non-iterative solver for wall-bounded domains. The PaScaL_TDMA library is used to solve batched tridiagonal systems, which are partitioned according to domain decomposition. For parallel fast Fourier transform in a Fourier Poisson solver for pressure, two transpose schemes with different communicator sizes are proposed and compared according to the number of cores. An explicit intermediate aggregation scheme for MPI-IO is suggested to reduce the number of processes that simultaneously take part in parallel IO. The overall implementation was evaluated using the Nurion cluster system of the Korea Institute of Science and Technology Information with detailed performance profiling. With the efficient communication, the results showed good strong and weak scalability up to 131,072 cores. The file IO performance improved with the suggested MPI-IO scheme using explicit aggregation, especially when many processes are involved in parallel IO. The solver was verified by conducting large-scale differentially heated vertical convection flow simulations under the Oberbeck–Boussinesq (OB) approximation. The capability of the solver to capture and quantitatively describe the non-OB effect was demonstrated via glycerol Rayleigh–Bénard convection flow simulations. Program summaryProgram Title: PaScaL_TCSCPC Library link to program files:https://doi.org/10.17632/3b2bysrcxm.1Developer's repository link:https://github.com/MPMC-LabLicensing provisions: MITProgramming language: Fortran90Nature of problem: Solving the incompressible Navier-Stokes equations for natural convective flow considering the non-Oberbeck–Boussinesq effect, which considers the temperature dependence of fluids properties, in three-dimensional channel domain.Solution method: PaScaL_TCS uses the monolithic projection-based method with staggered time discretization (MPM-STD) [1]. It employs the Crank-Nicolson scheme along with staggered time which evaluates scalar and vector variables at the (n+1/2) and (n+1) time level, respectively. PaScaL_TDMA [2] is used to solve multiple tridiagonal matrix systems for decoupled energy and momentum equations in a distributed memory system. The Poisson equation solver for the pressure-correction scheme is based on Fourier diagonalization. MPI_Alltoallw is used for parallel FFT. PaScaL_TCS supports a single-file-based parallel IO for efficient data generation in large-scale computing.

A conforming auxiliary space preconditioner for the mass conserving stress-yielding method.

We are studying the efficient solution of the system of linear equations stemming from the mass conserving stress-yielding (MCS) discretization of the Stokes equations. We perform static condensation to arrive at a system for the pressure and velocity unknowns. An auxiliary space preconditioner for the positive definite velocity block makes use of efficient and scalable solvers for conforming Finite Element spaces of low order and is analyzed with emphasis placed on robustness in the polynomial degree of the discretization. Numerical experiments demonstrate the potential of this approach and the efficiency of the implementation.

Parallel simulations for fast‐moving landslides: Space‐time mesh adaptation and sharp tracking of the wetting front

AbstractWe propose a highly scalable solver for a two‐dimensional depth‐integrated fluid dynamic model in order to simulate flow‐like landslides, such as debris or mud flows. The governing equations are discretized on quadtree meshes by means of a two‐step second‐order Taylor–Galerkin scheme, enriched by a suitable flux correction in order to avoid spurious oscillations, in particular near discontinuities and close to the wetting‐drying interface. A mesh adaptation procedure based on a gradient‐recovery a posteriori error estimator allows us to efficiently deal with a discretization of the domain customized to the phenomenon under investigation. Moreover, we resort to an adaptive scheme also in time to prevent filtering out the landslide dynamics, and to an interface tracking algorithm to avoid an excessive refinement in noninterfacial regions while preserving details along the wetting‐drying front. Finally, after verifying the performance of the proposed numerical method on idealized settings, we carry out a scalability analysis of the code both on idealized and real scenarios, to check the efficiency of the overall implementation.

Scalable computational algorithms for geospatial COVID-19 spread using high performance computing.

A nonlinear partial differential equation (PDE) based compartmental model of COVID-19 provides a continuous trace of infection over space and time. Finer resolutions in the spatial discretization, the inclusion of additional model compartments and model stratifications based on clinically relevant categories contribute to an increase in the number of unknowns to the order of millions. We adopt a parallel scalable solver that permits faster solutions for these high fidelity models. The solver combines domain decomposition and algebraic multigrid preconditioners at multiple levels to achieve the desired strong and weak scalabilities. As a numerical illustration of this general methodology, a five-compartment susceptible-exposed-infected-recovered-deceased (SEIRD) model of COVID-19 is used to demonstrate the scalability and effectiveness of the proposed solver for a large geographical domain (Southern Ontario). It is possible to predict the infections for a period of three months for a system size of 186 million (using 3200 processes) within 12 hours saving months of computational effort needed for the conventional solvers.

Arbitrary-Shape Dielectric Particles Interacting in the Linearized Poisson-Boltzmann Framework: An Analytical Treatment.

This work considers the interaction of two dielectric particles of arbitrary shape immersed into a solvent containing a dissociated salt and assuming that the linearized Poisson-Boltzmann equation holds. We establish a new general spherical re-expansion result which relies neither on the conventional condition that particle radii are small with respect to the characteristic separating distance between particles nor on any symmetry assumption. This is instrumental in calculating suitable expansion coefficients for the electrostatic potential inside and outside the objects and in constructing small-parameter asymptotic expansions for the potential, the total electrostatic energy, and forces in ascending order of Debye screening. This generalizes a recent result for the case of dielectric spheres to particles of arbitrary shape and builds for the first time a rigorous (exact at the Debye-Hückel level) analytical theory of electrostatic interactions of such particles at arbitrary distances. Numerical tests confirm that the proposed theory may also become especially useful in developing a new class of grid-free, fast, highly scalable solvers.

Variational satisfiability solving: efficiently solving lots of related SAT problems

Incremental satisfiability (SAT) solving is an extension of classic SAT solving that enables solving a set of related SAT problems by identifying and exploiting shared terms. However, using incremental solvers effectively is hard since performance is sensitive to the input order of subterms and results must be tracked manually. For analyses that generate sets of related SAT problems, such as those in software product lines, incremental solvers are either not used or their use is not clearly described in the literature. This paper translates the ordering problem to an encoding problem and automates the use of incremental solving. We introduce variational SAT solving, which differs from incremental solving by accepting all related problems as a single variational input and returning all results as a single variational output. Variational solving syntactically encodes differences in related SAT problems as local points of variation. With this syntax, our approach automates the interaction with the incremental solver and enables a method to automatically optimize sharing in the input. To evaluate these ideas, we formalize a variational SAT algorithm, construct a prototype variational solver, and perform an empirical analysis on two real-world datasets that applied incremental solvers to software evolution scenarios. We show, assuming a variational input, that the prototype solver scales better for these problems than four off-the-shelf incremental solvers while also automatically tracking individual results.

Robust and efficient primal-dual Newton-Krylov solvers for viscous-plastic sea-ice models

We present a Newton-Krylov solver for a viscous-plastic sea-ice model. This constitutive relation is commonly used in climate models to describe the material properties of sea ice. Due to the strong nonlinearity introduced by the material law in the momentum equation, the development of fast, robust and scalable solvers is still a substantial challenge. In this paper, we propose a novel primal-dual Newton linearization for the implicitly-in-time discretized momentum equation. Compared to existing methods, it converges faster and more robustly with respect to mesh refinement, and thus enables numerically converged sea-ice simulations at high resolutions. Combined with an algebraic multigrid-preconditioned Krylov method for the linearized systems, which contain strongly varying coefficients, the resulting solver scales well and can be used in parallel. We present experiments for two challenging test problems and study solver performance for problems with up to 8.4 million spatial unknowns.

Fully scalable solver for frequency-domain visco-elastic wave equations in 3D heterogeneous media: A controllability approach

We develop a controllability strategy for the computation of frequency-domain solutions of the 3D visco-elastic wave equation, in the perspective of seismic imaging applications. We generalize the controllability results for such equations beyond the sound-soft scattering (obstacle) problem. We detail the conjugate gradient implementation and show how an inner elliptic problem needs to be solved to compute the Riesz representative of the gradient at each iteration. We select a spectral-element spatial discretization and a fourth-order Runge-Kutta time discretization. We implement the controllability method in the framework of the SEM46 full waveform modeling and inversion software, to inherit for its excellent scalability which relies on an efficient domain decomposition algorithm. We perform a series of numerical experiments to validate the approach and illustrate its scalability up to more than fifteen hundred cores. In this case, with an elapsed time of less than 50 minutes, we solve a problem on a cubic domain containing up to 160 wavelengths in each direction, involving more than 1.7 billion unknowns.

Scalable multi-level deflation preconditioning for highly indefinite time-harmonic waves

Recent research efforts aimed at iteratively solving time-harmonic waves have focused on a broad range of techniques to accelerate convergence. In particular, for the famous Helmholtz equation, deflation techniques have been studied to accelerate the convergence of Krylov subspace methods. In this work, we extend the two-level deflation method to a multilevel deflation method for (heterogeneous) Helmholtz and elastic wave problems. By using higher-order deflation vectors, we show that up to the level where the coarse-grid linear systems remain indefinite, the near-zero eigenvalues of the these coarse-grid operators remain aligned with the near-zero eigenvalues of the fine-grid operator, keeping the spectrum of the preconditioned system away from the origin. Combining this with the well-known CSLP-preconditioner, we obtain a scalable solver for the highly indefinite linear systems. This can be attributed to a close to wave number independent convergence and an optimal use of the CSLP-preconditioner on the indefinite levels. There, we approximate the CSLP-preconditioner, while allowing the complex shift to be small, by using inner Bi-CGSTAB iterations instead of a multigrid F-cycle. The proposed method shows very promising results for the more challenging two- and three-dimensional heterogeneous time-harmonic wave problems.

Combining p-multigrid and Multigrid Reduction in Time methods to obtain a scalable solver for Isogeometric Analysis

The use of sequential time integration schemes becomes more and more the bottleneck within large-scale computations due to a stagnation of processor’s clock speeds. In this study, we combine the parallel-in-time Multigrid Reduction in Time method with a p-multigrid method to obtain a scalable solver specifically designed for Isogeometric Analysis. Numerical results obtained for two- and three-dimensional benchmark problems show the overall scalability of the proposed method on modern computer architectures and a significant improvement in terms of CPU timings compared to the use of standard spatial solvers.

GPU-accelerated scalable solver with bit permutated cyclic-min algorithm for quadratic unconstrained binary optimization

A wide range of combinatorial optimization problems can be reduced to the Ising model, and equivalently the quadratic unconstrained binary optimization (QUBO) problem. Thus, in recent years, researchers have proposed to solve QUBO on FPGAs, GPUs, and special-purpose processors. The adaptive bulk search (ABS) is a previously-proposed framework for solving QUBO in parallel on multiple GPUs. In the ABS, a CPU host performs a GA-based global search while GPUs asynchronously perform many local searches in parallel. An original ABS adopts a simple local search algorithm called a cyclic-min algorithm, which does not use pseudo random numbers. However, the lack of randomness may cause a potential drawback of restricted bit-flipping operations in a local search. To avoid this drawback, this paper proposes a cyclic-min algorithm with randomly-generated multiple bit permutations, which enables a more effective local search with random number generation in CPUs (not in GPUs). Furthermore, this paper introduces a scalable implementation of the ABS with MPI and OpenMP. Our experimental results on TSUBAME3.0 show that the solution quality improves and the throughput linearly increases as the number of GPUs increases; with 256 GPUs, it evaluates 20.1×1012 solutions per second.

Comparison of a finite volume and two Lattice Boltzmann solvers for swirled confined flows

A finite volume and two Lattice–Boltzmann unsteady, flow solvers using LES (Large Eddy Simulation) were compared in a swirling flow configuration, typical of aeronautical combustion chambers. Numerical results were validated against experimental data collected at EM2C laboratory by comparing pressure losses, mean and RMS velocity profiles on multiple planes and axial velocity spectra. Meshes and the overall numerical setups were individually adjusted for each code to obtain the targeted accuracy before comparing CPU efficiencies.Results confirm that the three LES codes provide high fidelity results, much better than usual RANS especially in terms of RMS data. The analysis of CPU performances shows that LBM (Lattice–Boltzmann Method) solvers are faster than the finite volume solver, even if CPU efficiencies remains of the same order of magnitude. In addition, strong scaling tests from 36 to 900 cores show that the finite volume solver scales more efficiently than the LBM codes, specially when the number of grid points per core is not sufficient.