Solver Library Research Articles

EllipticForest: A Direct Solver Library for Elliptic Partial Differential Equations on Adaptive Meshes

EllipticForest: A Direct Solver Library for Elliptic Partial Differential Equations on Adaptive Meshes

Hybrid programming-model strategies for GPU offloading of electronic structure calculation kernels.

To address the challenge of performance portability and facilitate the implementation of electronic structure solvers, we developed the basic matrix library (BML) and Parallel, Rapid O(N), and Graph-based Recursive Electronic Structure Solver (PROGRESS) library. The BML implements linear algebra operations necessary for electronic structure kernels using a unified user interface for various matrix formats (dense and sparse) and architectures (CPUs and GPUs). Focusing on density functional theory and tight-binding models, PROGRESS implements several solvers for computing the single-particle density matrix and relies on BML. In this paper, we describe the general strategies used for these implementations on various computer architectures, using OpenMP target functionalities on GPUs, in conjunction with third-party libraries to handle performance critical numerical kernels. We demonstrate the portability of this approach and its performance in benchmark problems.

Accurate error estimation for model reduction of nonlinear dynamical systems via data-enhanced error closure

Accurate error estimation is crucial in model order reduction, to obtain small reduced-order models as well as to certify their accuracy when deployed in downstream applications such as digital twins. In existing a posteriori error estimation approaches, knowledge about the time integration scheme is mandatory, e.g., the residual-based error estimators proposed for the reduced basis method. This poses a challenge when automatic ordinary differential equation solver libraries are used to perform the time integration. To address this, we present a data-enhanced approach for a posteriori error estimation. Our new formulation enables residual-based error estimators to be independent of any time integration method. To achieve this, we introduce a corrected reduced-order model which takes into account a data-driven closure term for improved accuracy. The closure term, subject to mild assumptions, is related to the local truncation error of the corresponding time integration scheme. We propose efficient computational schemes for approximating the closure term, at the cost of a modest amount of training data. Furthermore, the new error estimator is incorporated within a greedy process to obtain parametric reduced-order models. Numerical results on three different systems show the accuracy of the proposed error estimation approach and its ability to produce ROMs that generalize well.

HAZniCS – Software Components for Multiphysics Problems

We introduce the software toolbox HAZniCS for solving interface-coupled multiphysics problems. HAZniCS is a suite of modules that combines the well-known FEniCS framework for finite element discretization with solver and graph library HAZmath. The focus of this article is on the design and implementation of robust and efficient solver algorithms which tackle issues related to the complex interfacial coupling of the physical problems often encountered in applications in brain biomechanics. The robustness and efficiency of the numerical algorithms and methods is shown in several numerical examples, namely the Darcy-Stokes equations that model the flow of cerebrospinal fluid in the human brain and the mixed-dimensional model of electrodiffusion in the brain tissue.

Automated translation and accelerated solving of differential equations on multiple GPU platforms

We demonstrate a high-performance vendor-agnostic method for massively parallel solving of ensembles of ordinary differential equations (ODEs) and stochastic differential equations (SDEs) on GPUs. The method is integrated with a widely used differential equation solver library in a high-level language (Julia’s DifferentialEquations.jl) and enables GPU acceleration without requiring code changes by the user. Our approach achieves state-of-the-art performance compared to hand-optimized CUDA-C++ kernels while performing 20–100× faster than the vectorizing map (vmap) approach implemented in JAX and PyTorch. Performance evaluation on NVIDIA, AMD, Intel, and Apple GPUs demonstrates performance portability and vendor-agnosticism. We show composability with MPI to enable distributed multi-GPU workflows. The implemented solvers are fully featured – supporting event handling, automatic differentiation, and incorporation of datasets via the GPU’s texture memory – allowing scientists to take advantage of GPU acceleration on all major current architectures without changing their model code and without loss of performance. We distribute the software as an open-source library, DiffEqGPU.jl.

COMPARATIVE ANALYSIS OF SOFTWARE SYSTEMS FOR HYDRAULIC TURBINE FLOW SIMULATION

A comprehensive review of modern software complexes used for calculating spatial flow in hydraulic turbine flow parts was conducted. The widely used software system, Ansys, was analyzed. An overview of Ansys was provided, including its history, popularity within the CFD community, key features, and capabilities for analyzing the flow parts of hydraulic turbines. The preprocessing tools, solver parameters, post-processing functions, and visualization capabilities of Ansys were described. The advantages and limitations of using Ansys for calculating spatial flow in hydraulic turbine flow parts were analyzed. The open-source CFD software complex, OpenFOAM, was discussed. The main functions and capabilities of the OpenFOAM program were described. Information about solver libraries, meshing capabilities, advantages, and limitations for analyzing hydraulic turbines was presented, along with insights into the support from the scientific community and resources available to OpenFOAM users. SolidWorks FlowSimulation, which integrates with SolidWorks software, was examined. The unique features of SolidWorks FlowSimulation for analyzing spatial flow in hydraulic turbines were highlighted. The possibilities of CAD integration and the advantages of accurate geometric models were discussed. The capabilities of parametric analysis were explored, and the advantages and limitations of using SolidWorks FlowSimulation for calculating spatial flow in hydraulic turbine flow parts were analyzed. A comparison of the three software complexes was conducted based on their capabilities, ease of use, accuracy, computational resources required, and cost. An assessment of the advantages and disadvantages of each program was provided, along with recommendations for choosing the most suitable program based on specific use cases, objectives, and user requirements. This article serves as a valuable resource for engineers, designers, and researchers seeking insights into the available software systems for analyzing hydraulic turbine flow parts. It enables them to make informed decisions in selecting the most suitable software system based on their specific requirements, ultimately contributing to the optimization of hydraulic turbine performance and efficiency.

ARKODE: A Flexible IVP Solver Infrastructure for One-step Methods

We describe the ARKODE library of one-step time integration methods for ordinary differential equation (ODE) initial-value problems (IVPs). In addition to providing standard explicit and diagonally implicit Runge–Kutta methods, ARKODE supports one-step methods designed to treat additive splittings of the IVP, including implicit-explicit (ImEx) additive Runge–Kutta methods and multirate infinitesimal (MRI) methods. We present the role of ARKODE within the SUNDIALS suite of time integration and nonlinear solver libraries, the core ARKODE infrastructure for utilities common to large classes of one-step methods, as well as its use of “time stepper” modules enabling easy incorporation of novel algorithms into the library. Numerical results show example problems of increasing complexity, highlighting the algorithmic flexibility afforded through this infrastructure, and include a larger multiphysics application leveraging multiple algorithmic features from ARKODE and SUNDIALS.

Scalable Parallel Nonlinear Optimization with PyNumero and Parapint

We describe PyNumero, an open-source, object-oriented programming framework in Python that supports rapid development of performant parallel algorithms for structured nonlinear programming problems (NLP’s) using the Message Passing Interface (MPI). PyNumero provides three fundamental building blocks for developing NLP algorithms: a fast interface for calculating first and second derivatives with the AMPL Solver Library (ASL), a number of interfaces to efficient linear solvers, and block-structured vectors and matrices based on NumPy, SciPy, and MPI that support distributed parallel storage and computation. PyNumero’s design enables efficient, parallel algorithm development using high-level Python syntax while keeping expensive numerical calculations in fast, compiled implementations based on languages like C and Fortran. To demonstrate the utility of PyNumero, we also present Parapint, a Python package built on PyNumero for parallel solution of dynamic optimization problems. Parapint includes a parallel interior-point solver based on Schur-Complement decomposition. We illustrate the effectiveness of PyNumero for developing parallel algorithms with both code examples and scalability analyses for parallel matrix-vector dot products, parallel solution of structured systems of linear equations using Schur-Complement decomposition, and the parallel solution of a two-dimensional PDE optimal control problem. Our numerical results show nearly perfect scaling to more than 1,000 cores for large matrix-vector dot products and structured linear systems. Moreover, we obtain more than 354 times speedup for the optimal control example. History: Accepted by Alice Smith, EIC/Ted Ralphs, Area Editor/Software Tools. Funding: This work was funded in part by the Institute for the Design of Advanced Energy Systems (IDAES) with funding from the Office of Fossil Energy, Cross-Cutting Research, U.S. Department of Energy. This work was also funded by Sandia National Laboratories Laboratory Directed Research and Development (LDRD) program. Supplemental Material: The software that supports the findings of this study is available within the paper and its Supplemental Information ( https://pubsonline.informs.org/doi/suppl/10.1287/ijoc.2023.1272 ) as well as from the IJOC GitHub software repository ( https://github.com/INFORMSJoC/2021.0285 ) at ( http://dx.doi.org/10.5281/zenodo.7192328 ).

PANTR: A Proximal Algorithm With Trust-Region Updates for Nonconvex Constrained Optimization

This work presents PANTR, an efficient solver for nonconvex constrained optimization problems, that is well-suited as an inner solver for an augmented Lagrangian method. The proposed scheme combines forward-backward iterations with solutions to trust-region subproblems: the former ensures global convergence, whereas the latter enables fast update directions. We discuss how the algorithm is able to exploit exact Hessian information of the smooth objective term through a linear Newton approximation, while benefiting from the structure of box-constraints or ℓ1-regularization. An open-source C++ implementation of PANTR is made available as part of the NLP solver library ALPAQA. Finally, the effectiveness of the proposed method is demonstrated in nonlinear model predictive control applications.

Enabling New Flexibility in the SUNDIALS Suite of Nonlinear and Differential/Algebraic Equation Solvers

In recent years, the SUite of Nonlinear and DIfferential/ALgebraic equation Solvers (SUNDIALS) has been redesigned to better enable the use of application-specific and third-party algebraic solvers and data structures. Throughout this work, we have adhered to specific guiding principles that minimized the impact to current users while providing maximum flexibility for later evolution of solvers and data structures. The redesign was done through the addition of new linear and nonlinear solvers classes, enhancements to the vector class, and the creation of modern Fortran interfaces. The vast majority of this work has been performed “behind-the-scenes,” with minimal changes to the user interface and no reduction in solver capabilities or performance. These changes allow SUNDIALS users to more easily utilize external solver libraries and create highly customized solvers, enabling greater flexibility on extreme-scale, heterogeneous computational architectures.

SimELIT: A Novel GUI-Based Comprehensive Ion Trajectory Simulation Software for Mass Spectrometry.

Ion trajectory simulation in mass spectrometry systems from injection to detection is technically challenging but very important for better understanding the ion dynamics in instrument development. Here, we present SimELIT (Simulator of Eulerian and Lagrangian Ion Trajectories), a novel ion trajectory simulation platform. SimELIT is built upon a suite of multiphysics solvers compiled into OpenFOAM (an open-source numerical solver library particularly used for computational mechanics), with a simple web-based graphical user interface (GUI) allowing users to define the details of OpenFOAM cases and run simulations. SimELIT is a modular program and can provide extensions of physics (e.g., gas flows, electrodynamic fields) and thus enable ion trajectory simulations from the ion source to detector. The current version (SimELIT) provides two numerical solvers for ion trajectory simulations─(1) a Lagrangian particle tracker in vacuum and (2) a Eulerian ion density solver in background gas in the presence of electric fields. Here, we describe the architecture of SimELIT, including its use of Docker and the React Framework, and demonstrate the computation of ion trajectories of multiple m/z values in a static/linear voltage drop in vacuum (across a 1 m long flight tube). Further, the drift motion of ions under 1 Torr pressure conditions in a static background (N2) gas through a 20 V/cm static electric field is shown. The results produced from SimELIT were compared with SIMION and theoretical estimates. In addition, we report the computation of ion trajectories in electrodynamic fields within a planar FAIMS device operating at atmospheric pressure.

Toward a Lattice Boltzmann Method for Solids—Application to Static Equilibrium of Isotropic Materials

This work presents a novel method for simulating the behavior of solid objects with the Lattice Boltzmann Method (LBM). To introduce and validate our proposed framework, comparative studies are performed for computing the static equilibrium of isotropic materials. Remembering that the LBM has strong theoretical foundations in the Boltzmann equation; this latter is firstly adjusted to solid motions, through its Boltzmann-Vlasov special case. This is indeed the case when combined with a suitable mean-field external force term to set a reliable solid framework. Secondly, a library is built and plugged on the top of the well-known Parallel Lattice Boltzmann Solver (PaLaBoS) library. Numerical implementations based on the previous equation of motion for solids are led in a non-intrusive manner so as to present results with an easy and flawless reproducibility. A newly designed Lattice Boltzmann Method for Solids (LBMS) is exhibited through a few key algorithms, showing the overall operation plus the major improvements. Efficiency, robustness and accuracy of the proposed approach are illustrated and contrasted with a commercial Finite Element Analysis (FEA) software. The obtained results reveal considerable potential concerning static and further dynamic simulations involving solid constitutive laws within the LBM formalism.

Surface topography effects on dynamic behavior of water droplet over a micro-structured surface using an improved-VOF based lattice Boltzmann method

A volume of fluid-based lattice Boltzmann method (VOF-LBM) is developed, using the open-source parallel lattice Boltzmann solver library (Palabos) to study the wettability of solid surfaces and include the droplets’ dynamic behavior. This improves the accuracy and diminishes limitations on the density and viscosity ratios and the outcoming spurious velocities. An experimental study is performed, and the numerical results are verified through comparison with experimental data. Next, to design novel hydrophobic surfaces with pillar structures, the effect of different pillar patterns and topographies on the contact angle of water droplets and the droplet motion over an inclined surface are investigated. Hydrophobic micro-pillar surfaces are designed to enhance the ease of droplets movement on the slope. The effect of contact angle, wetted area, and surface energy on the sliding behavior of droplets is examined and presented. Results indicate that the new VOF-LBM method can simulate high density and viscosity ratios (1000 and 1.81×10-2 respectively), allowing real-world problems to be analyzed with appropriate accuracy. This reduces the spurious velocities produced at the fluids interface by 19 times compared to the Shan and Chen method. Results demonstrate that the contact angle of the water droplet on the surface with a contact angle of 115° can reach 160° by changing the surface topography. Also, a hydrophobic micro-pillar surface is designed in which the droplet moves nearly 12% faster than the flat surface.

Eulerian incompressible smoothed particle hydrodynamics on multiple GPUs

Recent advances in the development of Eulerian incompressible smoothed particle hydrodynamics (EISPH), such as high-order convergence and natural coupling with Lagrangian formulations, demonstrate its potential as a meshless alternative to traditional computational fluid dynamics (CFD) methods. This work aims to address one of the major outstanding limitations of EISPH, its relatively high computational cost, by providing an implementation that can be deployed on multiple graphics processing units (GPUs). To this end, a pre-existing multi-GPU version of the open-source Lagrangian weakly-compressible code DualSPHysics is converted to an EISPH formulation and integrated with an open-source multi-GPU multigrid solver (AmgX) to treat the pressure Poisson equation implicitly. The integration of AmgX within DualSPHysics presents a significant challenge, since AmgX is designed for distributed systems and therefore conflicts with the single-node shared memory design of the multi-GPU DualSPHysics code. The present implementation is validated against well-known test cases, showing excellent agreement with benchmark solutions and demonstrating second-order convergence. A detailed profiling and performance testing is also presented to investigate memory consumption and scaling characteristics. The results show approximately 87%–95% strong scaling efficiency and 92%–94% weak scaling efficiency in both 2D and 3D on up to four GPUs. Large spikes in memory consumption during the initialisation of the linear solver library are found to impede full utilisation of the device memory. Nevertheless, the present implementation is shown to permit problem sizes on the order of 69.3 million (2D) and 20.8 million (3D) particles on four GPUs (128 GB total device memory), which is beyond what has previously been reported for incompressible SPH on GPUs, where the PPE is treated implicitly, and demonstrates its potential as an alternative to traditional CFD methods.

Modelling and performance analysis of a single stage open absorption heat pump system in aspen plus using aqueous LiBr and HCOOK: A waste moist heat recovery application

Open absorption heat pump (OAHP) system is more viable option to recover waste heat from moist gas as compared to the traditional condensation methods. This promising technology has great potential for latent heat recovery from moist gas, drying process in paper and other industrial heating applications. This study presents the process modelling and comparative analysis of OAHP system in Aspen Plus using two different solutions by adopting part regeneration technique. The promising potassium formate-water (HCOOK/H2O) which has lower causticity, lower costs and better crystallization characteristics is used as an alternative to the caustic lithium bromide-water (LiBr/H2O) solution in this study. Process model of the system is established in Aspen Plus and, the properties validity is confirmed with published experimental and Engineering Equation Solver (EES) library data. A detailed comparative parametric study is carried out to evaluate the effect of influencing parameters on coefficient of performance (COP), water recovery (φ) and heat recovery (ζ) efficiencies. The performance of OAHP system is found to be very similar using different concentrations as 2.13 COP value for 50% LiBr/H2O and 2.19 for 70% HCOOK/H2O solution over design conditions. Similarly, φ is found to be 0.701, 0.688 while ζ as 0.716 and 0.705 for both the absorbents. Moreover, the system’s operational concentration range is 45-61.3% for LiBr/H2O and 55-82.1% for HCOOK/H2O at 135 °C regeneration heat input. Potassium formate solution having quite similar properties to the aqueous lithium bromide is also confirmed to have similar performance trends using 50% and 70% concentrations.

Effects of CO2 binary mixtures on pipeline performance for Carbon Capture and Storage

Captured CO2 from the manufacturing industries and power generation sources are not pure and may contain usually some amount of impure components, which affect the flow dynamics of the CO2 in the pipelines. Major component of impure components includes N2, H2, H2S, O2, Ar etc., depending on the capture technology used for the removal of CO2 streams. To design efficient CO2 pipeline transportation systems, it is imperative to understand the effect of these impure components on the flow behaviour. The simulation of the effect of CO2 binary mixture on pipeline performance is carried out in this study. The steady-state flow in pipeline is described by a set of parabolic mass, momentum and energy conservation equations. To solve the set of equations subject to the boundary and the inlet conditions of the pipeline, the non-linear algebraic solver library DNSQE in Fortran is used to study the behaviour of 90 vol% CO2, and 10 vol% single impure component (N2, H2, H2S, O2 and Ar). The phase diagram and pertinent fluid properties are calculated using the Peng-Robinson (PR) equation of state implemented in NIST Reference Fluid Thermodynamics and Transport Properties database program REFPROP. The results reveal that the bubble point curves were lifted up to higher pressures by impure components. The changes in the bubble point would affect the operating conditions of CO2 pipeline. The presence of impure components also changed the properties of the flowing fluid. From the results, H2 shows the most significant impact on the performance of CO2 pipeline as the highest pressure drop by ca. 0.18 bar/km which has smaller density and greater flow velocity against other CO2 streams. The results also show that Ar had the mildest effect followed by O2, N2 and CO. In the case of H2S, the flow parameters are found to be nearly identical to those for pure CO2. The pressure drop for binary mixture contains H2S shows the same value with pure CO2 stream by ca. 0.13 bar/km.

The Tailored CFD Package ‘containmentFOAM’ for Analysis of Containment Atmosphere Mixing, H2/CO Mitigation and Aerosol Transport

The severe reactor accident at Fukushima Daiichi Nuclear Power Plant (2011) has confirmed the need to understand the flow and transport processes of steam and combustible gases inside the containment and connected buildings. Over several years, Computational Fluid Dynamics (CFD) models, mostly based on proprietary solvers, have been developed to provide highly resolved insights; supporting the assessment of effectiveness of safety measures and possible combustion loads challenging the containment integrity. This paper summarizes the design and implementation of containmentFOAM, a tailored solver and model library based on OpenFOAM®. It is developed in support of Research & Development related to containment flows, mixing processes, pressurization, and assessment of passive safety systems. Based on preliminary separate-effect verification and validation results, an application oriented integral validation case is presented on the basis of an experiment on gas mixing and H2 mitigation by means of passive auto-catalytic recombiners in the THAI facility (Becker Technologies, Eschborn, Germany). The simulation results compare well with the experimental data and demonstrate the general applicability of containmentFOAM for technical scale analysis. Concluding the paper, the strategy for dissemination of the code and measures implemented to minimize potential user errors are outlined.

Effects of CO2 binary mixtures on pipeline performance for Carbon Capture and Storage

Captured CO2 from the manufacturing industries and power generation sources are not pure and may contain usually some amount of impure components, which affect the flow dynamics of the CO2 in the pipelines. Major component of impurities includes N2, H2, H2S, O2, Ar etc., depending on the capture technology used for the removal of CO2 [1]. To design efficient CO2 pipeline transportation systems, it is imperative to understand the effect of these impure components on the flow behaviour. The simulation of the effect of CO2 binary mixture on pipeline performance is carried out in this study. The steady-state flow in pipeline is described by a set of parabolic mass, momentum and energy conservation equations. To solve the set of equations subject to the boundary and the inlet conditions of the pipeline, the nonlinear algebraic solver library DNSQE in Fortran is used to study the behaviour of 90 vol % CO2, and 10 vol % single impure component (N2, H2, H2S, O2 and Ar). The pertinent fluid properties are calculated using the Peng-Robinson (PR) equation of state. The presence of impure components changed the properties of the flowing fluid. From the results, H2 shows the most significant impact on the performance of CO2 pipeline as the highest pressure drop has smaller density and greater flow velocity against other CO2 streams. The results also show that Ar had the mildest effect followed by O2, N2 and CO. In the case of H2S, the flow parameters are found to be nearly identical to those for pure CO2.

Phase-field modeling of multivariant martensitic transformation at finite-strain: Computational aspects and large-scale finite-element simulations

Large-scale 3D martensitic microstructure evolution problems are studied using a finite-element discretization of a finite-strain phase-field model. The model admits an arbitrary crystallography of transformation and arbitrary elastic anisotropy of the phases, and incorporates Hencky-type elasticity, a penalty-regularized double-obstacle potential, and viscous dissipation. The finite-element discretization of the model is performed in Firedrake and relies on the PETSc solver library. The large systems of linear equations arising are efficiently solved using GMRES and a geometric multigrid preconditioner with a carefully chosen relaxation. The modeling capabilities are illustrated through a 3D simulation of the microstructure evolution in a pseudoelastic CuAlNi single crystal during nano-indentation, with all six orthorhombic martensite variants taken into account. Robustness and a good parallel scaling performance have been demonstrated, with the problem size reaching 150 million degrees of freedom.

FLUPS: A Fourier-Based Library of Unbounded Poisson Solvers

A Fourier-based Library of Unbounded Poisson Solvers (FLUPS) for 2D and 3D homogeneous distributed grids is presented. It is designed to handle every possible combination of periodic, symmetric, semi-unbounded and fully unbounded boundary conditions for the Poisson equation on rectangular domains with uniform resolution. FLUPS leverages a dedicated implementation of 3D Fourier transforms to solve the Poisson equation using Green's functions, in a fast and memory-efficient way. Several Green's functions are available, optionally with explicit regularization, spectral truncation, or using lattice Green's functions, and provide verified convergence orders from 2 to spectral-like. The algorithm depends on the FFTW library to perform 1D transforms, while Message Passing Interface (MPI) communications enable the required remapping of data in memory. For the latter operation, a first available implementation resorts to the standard all-to-all routines. A second implementation, featuring non-blocking and persistent point-to-point communications, is however shown to be more efficient in a majority of cases and especially while taking advantage of the shared memory parallelism with OpenMP. The scalability of the algorithm, aimed at massively parallel architectures, is demonstrated up to 73 720 cores. The results obtained with three different supercomputers show that the weak efficiency remains above 40\% and the strong efficiency above 30% when the number of cores is multiplied by 16, for typical problems. These figures are slightly better than those expected from a third party 3D Fast Fourier Transform (FFT) tool, with which a 20% longer execution time was also measured on average.