Solver Library Research Articles

PaScaL_TDMA: A library of parallel and scalable solvers for massive tridiagonal systems

The aim of this study is to devise an efficient and scalable computational procedure to solve the many tridiagonal systems in multi-dimensional partial differential equations. The modified Thomas algorithm and a newly designed communication scheme were used to reduce the communication overhead encountered while solving the many tridiagonal systems. Benchmark test results reveal an advantage of the proposed procedures compared to global all-to-all communication methods — a significantly reduced communication time that becomes more prominent for larger problem sizes and greater number of cores. The proposed computational procedures are fully implemented in an open-source library called Parallel and Scalable Library for TDMA (PaScaL_TDMA). Considering a three-dimensional heat conduction problem as a practical example, we obtain good strong and weak scalability results up to 262,144 computing cores on the KISTI Nurion cluster system, which, to the best of our knowledge, is the largest parallel simulation for solving tridiagonal systems. The potential of this library for large-scale substantive problems in physics is also demonstrated through direct numerical simulations of the Rayleigh–Bénard convection problem, which yielded excellent scalability and accurate results. Program summaryProgram Title: PaScaL_TDMACPC Library link to program files:http://dx.doi.org/10.17632/49z6fh94z3.1Developer’s repository link:https://github.com/MPMC-LabLicensing provisions: MITProgramming language: Fortran90Nature of problem: Tridiagonal systems for solving multi-dimensional partial differential equations.Solution method: The divide-and-conquer method is used to solve partitioned tridiagonal systems of equations in the distributed memory system. The partitioned tridiagonal systems of equations are transformed into modified sub-matrices using the modified Thomas algorithm [1]. Reduced tridiagonal systems are constructed by collecting the first and last rows of the modified submatrices from each computing core. The newly designed communication scheme based on MPI_Alltoallw accelerates to collect rows and construct the reduced tridiagonal systems. The solutions of the reduced tridiagonal systems are obtained via the sequential Thomas algorithm. Thereafter, the remaining unknowns in the modified sub-matrices are solved using solutions of the reduced tridiagonal systems.

A Feature-complete SPIKE Dense Banded Solver

This article presents a parallel, effective, and feature-complete recursive SPIKE algorithm that achieves near feature-parity with the standard linear algebra package banded linear system solver. First, we present a flexible parallel implementation of the recursive SPIKE scheme that aims at removing its original limitation that the number of cores/processors be restricted to powers of two. A new transpose solve option for SPIKE is then developed to satisfy a standard requirement of most numerical solver libraries. Finally, a pivoting recursive SPIKE strategy is presented as an alternative to the non-pivoting scheme to improve numerical stability. All these new enhancements lead to the release of a new black-box feature-complete SPIKE-OpenMP package that significantly improves upon the performance and scalability obtained with other state-of-the-art banded solvers.

ELSI — An open infrastructure for electronic structure solvers

Routine applications of electronic structure theory to molecules and periodic systems need to compute the electron density from given Hamiltonian and, in case of non-orthogonal basis sets, overlap matrices. System sizes can range from few to thousands or, in some examples, millions of atoms. Different discretization schemes (basis sets) and different system geometries (finite non-periodic vs. infinite periodic boundary conditions) yield matrices with different structures. The ELectronic Structure Infrastructure (ELSI) project provides an open-source software interface to facilitate the implementation and optimal use of high-performance solver libraries covering cubic scaling eigensolvers, linear scaling density-matrix-based algorithms, and other reduced scaling methods in between. In this paper, we present recent improvements and developments inside ELSI, mainly covering (1) new solvers connected to the interface, (2) matrix layout and communication adapted for parallel calculations of periodic and/or spin-polarized systems, (3) routines for density matrix extrapolation in geometry optimization and molecular dynamics calculations, and (4) general utilities such as parallel matrix I/O and JSON output. The ELSI interface has been integrated into four electronic structure code projects (DFTB+, DGDFT, FHI-aims, SIESTA), allowing us to rigorously benchmark the performance of the solvers on an equal footing. Based on results of a systematic set of large-scale benchmarks performed with Kohn–Sham density-functional theory and density-functional tight-binding theory, we identify factors that strongly affect the efficiency of the solvers, and propose a decision layer that assists with the solver selection process. Finally, we describe a reverse communication interface encoding matrix-free iterative solver strategies that are amenable, e.g., for use with planewave basis sets. Program summaryProgram title: ELSI InterfaceCPC Library link to program files:http://dx.doi.org/10.17632/473mbbznrs.1Licensing provisions: BSD 3-clauseProgramming language: Fortran 2003, with interface to C/C++External routines/libraries: BLACS, BLAS, BSEPACK (optional), EigenExa (optional), ELPA, FortJSON, LAPACK, libOMM, MPI, MAGMA (optional), MUMPS (optional), NTPoly, ParMETIS (optional), PETSc (optional), PEXSI, PT-SCOTCH (optional), ScaLAPACK, SLEPc (optional), SuperLU_DISTNature of problem: Solving the electronic structure from given Hamiltonian and overlap matrices in electronic structure calculations.Solution method: ELSI provides a unified software interface to facilitate the use of various electronic structure solvers including cubic scaling dense eigensolvers, linear scaling density matrix methods, and other approaches.

Siesta: Recent developments and applications.

A review of the present status, recent enhancements, and applicability of the Siesta program is presented. Since its debut in the mid-1990s, Siesta's flexibility, efficiency, and free distribution have given advanced materials simulation capabilities to many groups worldwide. The core methodological scheme of Siesta combines finite-support pseudo-atomic orbitals as basis sets, norm-conserving pseudopotentials, and a real-space grid for the representation of charge density and potentials and the computation of their associated matrix elements. Here, we describe the more recent implementations on top of that core scheme, which include full spin-orbit interaction, non-repeated and multiple-contact ballistic electron transport, density functional theory (DFT)+U and hybrid functionals, time-dependent DFT, novel reduced-scaling solvers, density-functional perturbation theory, efficient van der Waals non-local density functionals, and enhanced molecular-dynamics options. In addition, a substantial effort has been made in enhancing interoperability and interfacing with other codes and utilities, such as wannier90 and the second-principles modeling it can be used for, an AiiDA plugin for workflow automatization, interface to Lua for steering Siesta runs, and various post-processing utilities. Siesta has also been engaged in the Electronic Structure Library effort from its inception, which has allowed the sharing of various low-level libraries, as well as data standards and support for them, particularly the PSeudopotential Markup Language definition and library for transferable pseudopotentials, and the interface to the ELectronic Structure Infrastructure library of solvers. Code sharing is made easier by the new open-source licensing model of the program. This review also presents examples of application of the capabilities of the code, as well as a view of on-going and future developments.

Design and Verification of Tuned Liquid Column Dampers for High-rise Buildings Using CFD

BMT Fluid Mechanics conducted a comprehensive flow study on the effectiveness of a pair of Tuned Liquid Column Dampers (TLCDs) aimed at mitigating wind-induced accelerations in a high-rise building. The study successfully combined Computational Fluid Dynamics (CFD) with physical model tests performed at the 6 degree-of-freedom shake table facility located at the Department of Civil Engineering of the University of Bristol. The original TLCD concept discussed in this work was proposed at the early design stages of a slender 42-storey high-end residential building located in the Middle East. On-site full-scale measurements and High Frequency Force Balance (HFFB) wind tunnel tests performed by BMT Fluid Mechanics showed that the highest occupied floors could experience excessive wind-induced motion. Depending on the inherent damping of the finished structure, this motion had the potential to exceed standard occupant comfort criteria. A new CFD methodology was developed to assess the performance of the TLCD design using the multi-phase Volume-of-Fluid (VOF) solver and dynamic mesh motion libraries available in the CFD software tool HELYX®. The damping effects introduced by the array of internal porous baffles in the TLCD were approximated using a Darcy-Forchheimer porosity model. Comparisons of free-decay damping performance between CFD results and shake table experiments for a 1:20 scale model of a 3 baffles’ TLCD configuration were found to be satisfactory for design purposes. The results for the amplitude of the net forces acting on the TLCD, as well as the frequency response measured using the free-decay logarithmic decrement approach, confirmed that the proposed CFD methodology provides an accurate representation of real operating conditions. The same CFD approach was successfully applied to model the full-scale TLCD device. With the introduction of a pair of identical purposely designed TLCDs, a 30% increase performance of the structural response of the 42-storey building to wind loading excitation was achieved.

Prediction of the flow field and convective heating during space capsule reentry using an open source solver

A numerical investigation is carried out using an Open source solver for a hypersonic flow around a space capsule. Calculations are performed to test the reliability of the solver to simulate Martian reentry flows. The kinetic model and viscosity coefficients have been added to the solver library. A mesh study for the fore-body concluded that a medium mesh was sufficient to capture the shock standoff distance with remarkable accuracy, which leads us to believe that the numerical scheme used is quite robust. The study then focused on the effect of wall condition on the flow topology and proved once again that the solver could accurately take the influence of boundary conditions on flow parameters. Following the obtained results, hy2Foam can be considered as a reliable tool for numerical simulations of space capsules reentry.

Optimizations of the eigensolvers in the ELPA library

The solution of (generalized) eigenvalue problems for symmetric or Hermitian matrices is a common subtask of many numerical calculations in electronic structure theory or materials science. Depending on the scientific problem, solving the eigenvalue problem can easily amount to a sizeable fraction of the whole numerical calculation, and quite often is even the dominant part by far. For researchers in the field of computational materials science, an efficient and scalable solution of the eigenvalue problem is thus of major importance. The ELPA-library (Eigenvalue SoLvers for Petaflop-Applications) is a well-established dense direct eigenvalue solver library, which has proven to be very efficient and scalable up to very large core counts. It is in a wide-spread production use on a large variety of HPC systems worldwide, and is applied by many codes in the field of materials science. In this paper, we describe the latest optimizations of the ELPA-library for new HPC architectures of the Intel Skylake processor family with an AVX-512 SIMD instruction set, or for HPC systems accelerated with recent GPUs. Apart from those direct hardware-related optimizations, we also describe a complete redesign of the API in a modern modular way, which, apart from a much simpler and more flexible usability, leads to a new path to access system-specific performance optimizations. In order to ensure optimal performance for a particular scientific setting or a specific HPC system, the new API allows the user to influence in a straightforward way the internal details of the algorithms and of performance-critical parameters used in the ELPA-library. On top of that, we introduce an autotuning functionality, which allows for finding the best settings in a self-contained automated way, without the need of much user effort. In situations where many eigenvalue problems with similar settings have to be solved consecutively, the autotuning process of the ELPA-library can be done “on-the-fly”, without the need of preceding the simulation with an “artificial” autotuning step. Practical applications from materials science which rely on reaching a numerical convergence limit by so-called self-consistency iterations can profit from the on-the-fly autotuning. On some examples of scientific interest, simulated with the FHI-aims application, the advantages of the latest optimizations of the ELPA-library are demonstrated.

Pre-conditioning strategies to accelerate the convergence of iterative methods in multiphase flow simulations

A computational bottleneck during the solution to multiphase formulations of the incompressible Navier–Stokes equations is often during the implicit solution of the pressure-correction equation that results from operator-splitting methods. Since density is a coefficient in the pressure-correction equation, large variations or discontinuities among the phase densities greatly increase the condition number of the pressure-correction matrix and retard the convergence of iterative methods employed in its solution. To alleviate this shortcoming, the open-source multiphase code MFiX is interfaced with the linear solver library PETSc. Through an appropriate mapping of matrix and vector data structures between the two software packages, an access to a suite of robust, scalable, solver options in PETSc is obtained. Verification of the implementation is demonstrated through predictions that are identical to those obtained from MFiX’s native solvers for a class of single-phase and multiphase flow problems. For a low Reynolds number, flow over a cylinder case, applying Right Side Block Jacobi Preconditioning to the BiCGSTAB iterative solver in PETSc was faster than MFiX’s native solver. This speed-up increased with higher mesh resolution and for higher-order spatial discretizations. In a fluidized bed simulation, this solver–preconditioner combination resulted in a 25% decrease in solve time compared to MFiX’s native solver.

Description and evaluation of the Community Ice Sheet Model (CISM) v2.1

Abstract. We describe and evaluate version 2.1 of the Community Ice Sheet Model (CISM). CISM is a parallel, 3-D thermomechanical model, written mainly in Fortran, that solves equations for the momentum balance and the thickness and temperature evolution of ice sheets. CISM's velocity solver incorporates a hierarchy of Stokes flow approximations, including shallow-shelf, depth-integrated higher order, and 3-D higher order. CISM also includes a suite of test cases, links to third-party solver libraries, and parameterizations of physical processes such as basal sliding, iceberg calving, and sub-ice-shelf melting. The model has been verified for standard test problems, including the Ice Sheet Model Intercomparison Project for Higher-Order Models (ISMIP-HOM) experiments, and has participated in the initMIP-Greenland initialization experiment. In multimillennial simulations with modern climate forcing on a 4 km grid, CISM reaches a steady state that is broadly consistent with observed flow patterns of the Greenland ice sheet. CISM has been integrated into version 2.0 of the Community Earth System Model, where it is being used for Greenland simulations under past, present, and future climates. The code is open-source with extensive documentation and remains under active development.

Sundials/ML: Connecting OCaml to the Sundials Numeric Solvers

This paper describes the design and implementation of a comprehensive OCaml interface to the Sundials library of numeric solvers for ordinary differential equations, differential algebraic equations, and non-linear equations. The interface provides a convenient and memory-safe alternative to using Sundials directly from C and facilitates application development by integrating with higher-level language features, like garbage-collected memory management, algebraic data types, and exceptions. Our benchmark results suggest that the interface overhead is acceptable: the standard examples are rarely twice as slow in OCaml than in C, and often less than 50% slower. The challenges in interfacing with Sundials are to efficiently and safely share data structures between OCaml and C, to support multiple implementations of vector operations and linear solvers through a common interface, and to manage calls and error signalling to and from OCaml. We explain how we overcame these difficulties using a combination of standard techniques such as phantom types and polymorphic variants, and carefully crafted data representations.

Incremental SAT Library Integration Using Abstract Stobjs

We describe an effort to soundly use off-the-shelf incremental SAT solvers within ACL2 by modeling the behavior of a SAT solver library as an abstract stobj. The interface allows ACL2 programs to use incremental SAT solvers, and the abstract stobj model allows us to reason about the behavior of an incremental SAT library so as to show that algorithms implemented using it are correct, as long as the library is bug-free.

Supercritical steam outflow through a thin nozzle: forming a hollow jet

The process of hollow jet formation during steam outflow through a thin nozzle was studied; the water steam was initially at high pressure and in the supercritical state. The numerical solution for this problem was obtained with the sonicFoam solver library from the open-source CFD software package OpenFOAM in 2D axisymmetric formula-tion. The reliability of results is estimated by comparing two approaches for simulation of the dynamics of unloading wave propagating in the high-pressure nozzle using the OpenFOAM package with Peng−Robinson equation of state and numerical solution of a similar problem by method of through computation in the case of 1D planar approximation for perfect gas equation of state.

Parallel implementation of implicit finite element model with cohesive zones and collision response using CUDA

SummaryThe aim of this work is to efficiently implement the Park‐Paulino‐Roesler cohesive zone model with the objective of creating realistic high‐resolution simulations of material deformation, fracture, and postfracture behavior. Intrinsically, unstructured meshes can create more realistic fracture patterns in bulk material than structured meshes. Implicit methods, stable for much larger time steps, have greater potential to model both fracture and postfracture behavior without sacrificing speed of execution. Several technical contributions are presented, including (i) GPU‐accelerated implementation of the Park‐Paulino‐Roesler cohesive zone model, (ii) efficient creation of sparse matrix structure, and (iii) comparison of different unloading/reloading relations when using an implicit scheme. A potential‐based collision response scheme was implemented that allows one to model the interaction of fragmented material. Several test simulations are carried out to demonstrate the flexibility of the model and its ability to reproduce different materials under various loading conditions. Benchmarking results show that most of the computational time is spent by the third‐party solver library, meaning that other operations do not require optimization. The library is made available as open source.

DL_MG: A Parallel Multigrid Poisson and Poisson-Boltzmann Solver for Electronic Structure Calculations in Vacuum and Solution.

The solution of the Poisson equation is a crucial step in electronic structure calculations, yielding the electrostatic potential-a key component of the quantum mechanical Hamiltonian. In recent decades, theoretical advances and increases in computer performance have made it possible to simulate the electronic structure of extended systems in complex environments. This requires the solution of more complicated variants of the Poisson equation, featuring nonhomogeneous dielectric permittivities, ionic concentrations with nonlinear dependencies, and diverse boundary conditions. The analytic solutions generally used to solve the Poisson equation in vacuum (or with homogeneous permittivity) are not applicable in these circumstances, and numerical methods must be used. In this work, we present DL_MG, a flexible, scalable, and accurate solver library, developed specifically to tackle the challenges of solving the Poisson equation in modern large-scale electronic structure calculations on parallel computers. Our solver is based on the multigrid approach and uses an iterative high-order defect correction method to improve the accuracy of solutions. Using two chemically relevant model systems, we tested the accuracy and computational performance of DL_MG when solving the generalized Poisson and Poisson-Boltzmann equations, demonstrating excellent agreement with analytic solutions and efficient scaling to ∼109 unknowns and 100s of CPU cores. We also applied DL_MG in actual large-scale electronic structure calculations, using the ONETEP linear-scaling electronic structure package to study a 2615 atom protein-ligand complex with routinely available computational resources. In these calculations, the overall execution time with DL_MG was not significantly greater than the time required for calculations using a conventional FFT-based solver.

A parallel framework for a multipoint flux mixed finite element equation of state compositional flow simulator

Mathematical models of physical problems are becoming increasingly complex and computationally intensive. At the same time, computing hardware is becoming more parallelized with an increasing number of cores promoting simultaneous tasks. In this work, we present a parallel, equation of state (EOS), compositional flow simulator for evaluating CO2 sequestration, enhanced oil recovery techniques such as gas flooding, and other subsurface porous media applications. Using the multipoint flux mixed finite element (MFMFE) method for spatial discretization, it can handle complex reservoir geometries using general distorted hexahedral grid elements, as well as satisfy local mass conservation and compute accurate phase fluxes. A parallel framework for the MFMFE is presented that has been extended to the highly non-linear, EOS, compositional flow model. Much of the non-linearity is due to the local flash and stability calculations associated with interphase mass transfer and phase behavior. Parallel multigrid linear solver libraries such as HYPRE are utilized to solve the algebraic problems on each Newton step. We perform a variety of strong and weak parallel scaling studies up to 10 million elements and 1024 processors, and discuss possible load balancing issues.

A computationally efficient P1 radiation model for modern combustion systems utilizing pre-conditioned conjugate gradient methods

The iterative convergence of the P1 radiation model can be slow in optically thin scenarios when employing classical iterative methods. In order to remedy this shortcoming, an in-house P1 radiation model was interfaced with high performance, scalable, linear solver libraries. Next, the accuracies of P1 radiation model calculations was assessed by comparing its predictions against discrete ordinates (DO) model calculations for prototypical problems representative of modern combustion systems. Corresponding benchmark results were also included for comparison.Utilizing Pre-Conditioners (PC) to the Conjugate Gradients (CG) method, the convergence time of the P1 radiation model reduced by a factor of 30 for modest problem sizes and a factor of 70 for larger sized problems when compared against classical Gauss Seidel sweeps. Further, PC provided 50% computational savings compared to employing CG in a standalone mode. The P1 model calculation times were about 25–30% of the DO model calculation time. The time to solution also scaled linearly with an increase in problem size. The weighted sum of gray gases model employed in this study in conjunction with the P1 model provided good agreement against benchmark data with L2 error norms (defined relative to corresponding DO calculations) improving when isotropic intensities were prevalent.

Parameter study and CFD analysis of head on collision and dynamic behavior of two colliding ferrofluid droplets

In the present study a computational fluid dynamics approach is implemented to investigate the dynamic behavior of two freely suspended ferrofluid droplets under the effect of uniform magnetic field. The colliding droplets are tracked from their initial state to a new equilibrium state which will be obtained for the new produced droplet after the coalescence. During collision time the shape of each droplet and also the variations of their energies are investigated and discussed. For the simulation, a finite volume based solver is modified based on the open source solver library (OpenFOAM®) which is capable of coupling the flow field and magnetostatic equations. A coupled scheme of volume of fluid and level-set methods is applied for interface capturing. The simulation results are validated based on the available numerical and experimental results reported in the literature. In order to evaluate the effect of parameters on the magnetic induced coalescence of ferrofluid droplets a dimensional analysis based on the Buckingham π theorem is implemented. Based on the simulation results and dimensional analysis a new correlation is proposed which is capable to predict the collision time for two ferrofluid droplets in a relatively wide range of properties and operating conditions.

GR-AMRVAC code applications: accretion onto compact objects, boson stars versus black holes

In the close vicinity of a compact object strong gravity imprints its signature onto matter. Systems that contain at least one compact object are observed to exhibit extreme physical properties and typically emit highly energetic radiation. The nature of the compact objects that produce the strongest gravitational fields is to date not settled. General relativistic numerical simulations of fluid dynamics around black holes, neutron stars, and other compact objects such as boson stars (BSs) may give invaluable insights into this fundamental question. In order to study the behavior of fluid in the strong gravity regime of an arbitrary compact object we develop a new general relativistic hydrodynamics code. To this end we extend the existing versatile adaptive mesh refinement code MPI-AMRVAC into a general relativistic hydrodynamics framework and adapt it for the use of numerically given spacetime metrics. In the present article we study accretion flows in the vicinity of various types of BSs whose numerical metrics are calculated by the KADATH spectral solver library. We design specific tests to check the reliability of any code intending to study BSs and compare the solutions with those obtained in the context of Schwarzschild black holes. We perform the first ever general relativistic hydrodynamical simulations of gas accretion by a BS. The behavior of matter at small distances from the center of a BS differs notably from the black hole case. In particular we demonstrate that in the context of Bondi spherical accretion the mass accretion rate onto non-rotating BSs remains constant whereas it increases for Schwarzschild black holes. We also address the scenario of non-spherical accretion onto BSs and show that this may trigger mass ejection from the interior of the BS. This striking feature opens the door to forthcoming investigations regarding accretion-ejection flows around such types of compact objects.

Strongly coupled partitioned six degree-of-freedom rigid body motion solver with Aitken's dynamic under-relaxation

An implicit method of solving the six degree-of-freedom rigid body motion equations based on the second order Adams-Bashforth-Moulten method was utilised as an improvement over the leapfrog scheme by making modifications to the rigid body motion solver libraries directly. The implementation will depend on predictor-corrector steps still residing within the hybrid Pressure Implicit with Splitting of Operators - Semi-Implicit Method for Pressure Linked Equations (PIMPLE) outer corrector loops to ensure strong coupling between fluid and motion. Aitken's under-relaxation is also introduced in this study to optimise the convergence rate and stability of the coupled solver. The resulting coupled solver ran on a free floating object tutorial test case when converged matches the original solver. It further allows a varying 70%–80% reduction in simulation times compared using a fixed under-relaxation to achieve the required stability.

Laplace-domain waveform modeling and inversion for the 3D acoustic–elastic coupled media

Laplace-domain waveform inversion reconstructs long-wavelength subsurface models by using the zero-frequency component of damped seismic signals. Despite the computational advantages of Laplace-domain waveform inversion over conventional frequency-domain waveform inversion, an acoustic assumption and an iterative matrix solver have been used to invert 3D marine datasets to mitigate the intensive computing cost. In this study, we develop a Laplace-domain waveform modeling and inversion algorithm for 3D acoustic–elastic coupled media by using a parallel sparse direct solver library (MUltifrontal Massively Parallel Solver, MUMPS). We precisely simulate a real marine environment by coupling the 3D acoustic and elastic wave equations with the proper boundary condition at the fluid-solid interface. In addition, we can extract the elastic properties of the Earth below the sea bottom from the recorded acoustic pressure datasets. As a matrix solver, the parallel sparse direct solver is used to factorize the non-symmetric impedance matrix in a distributed memory architecture and rapidly solve the wave field for a number of shots by using the lower and upper matrix factors. Using both synthetic datasets and real datasets obtained by a 3D wide azimuth survey, the long-wavelength component of the P-wave and S-wave velocity models is reconstructed and the proposed modeling and inversion algorithm are verified. A cluster of 80CPU cores is used for this study.