Algebraic Solver Research Articles

The Aqueous Aluminium-Ion Battery: Optimising the Electrode Compression Ratio through Image-Based Modelling

The aqueous Al-ion battery potentially has improved safety and environmental advantages over the incumbent lithium-ion technology [1,2]. Using 3D carbon-felt matrix electrodes, the performance of these batteries can be optimised. Typically, the electrodes are compressed to improve electrical conductivity by improving contact between adjacent fibres, which also closes the electrolyte pores, reducing the ionic diffusivity [3]. In this work, we set out to understand the role of electrode compression on the interplay between these key performance parameters for the aq. Al-ion battery.Synchrotron X-ray computed tomography was used to examine the 3D morphology of the carbon felt electrodes under 11 different compression ratios [4]. A bespoke in-situ tensile/compression rig was used to measure displacement and loading for three electrode types: raw carbon felt, positive electrodes loaded with a copper hexacyanoferrate ink and negative electrodes loaded with TiO2 anatase powder. Data was acquired using two voxel sizes (330nm, and 540nm) to compare the effects of spatial resolution and imaged volume size on subsequent image analysis. The tomograms were reconstructed using Paganin phase retrieval, to improve the contrast of the weakly attenuating carbon, and a filtered back projection algorithm. A U-net convolutional neural network was trained on uncompressed, fully compressed and partially compressed (three volumes) data, and then used to fully segment all tomograms.A high-throughput, image-based model was then used to analyse how compression affects the porosity, tortuosity, volume-specific area, ionic diffusivity, and electrical conduction of the electrodes. A finite differences-based model was used to solve the equilibrium partial differential equations directly on the voxel datasets, with no additional regularisation [5]. The heterogeneity of the electrode samples is quantified by comparing the effect of representative elementary volume on the value of the computed parameters [6]. Finally, aq. Al-ion batteries are manufactured using the predicted optimal compression ratio and compared to the simulated results.This is the first in-situ study of the compression effect on the aqueous aluminium-ion battery, and the largest XCT-based modelling study known to the authors (99 XCT tomograms). This work aims to improve the understanding of the effect of manufacturing parameters on the aqueous Al-ion battery and other similar batteries, and ultimately, its performance.[1] Melzack, N. "Advancing battery design based on environmental impacts using an aqueous Al-ion cell as a case study." Scientific Reports 12, no. 1 (2022): 8911.[2] Melzack, Nicole, R. G. A. Wills, and Andrew Cruden. "Cleaner energy storage: Cradle-to-gate life cycle assessment of aluminum-ion batteries with an aqueous electrolyte." Frontiers in Energy Research 9 (2021): 699919.[3] Mckerracher, R. D., A. Holland, A. Cruden, and R. G. A. Wills. "Comparison of carbon materials as cathodes for the aluminium-ion battery." Carbon 144 (2019): 333-341.[4] Reinhard, Christina, Michael Drakopoulos, Sharif I. Ahmed, Hans Deyhle, Andrew James, Christopher M. Charlesworth, Martin Burt et al. "Beamline K11 DIAD: A new instrument for dual imaging and diffraction at Diamond Light Source." Journal of Synchrotron Radiation 28, no. 6 (2021): 1985-1995.[5] Le Houx, James, and Denis Kramer. "Openimpala: Open source image based parallisable linear algebra solver." SoftwareX 15 (2021): 100729.[6] Le Houx, James, Siul Ruiz, Daniel McKay Fletcher, Sharif Ahmed, and Tiina Roose. "Statistical Effective Diffusivity Estimation in Porous Media Using an Integrated On-site Imaging Workflow for Synchrotron Users." Transport in Porous Media 150, no. 1 (2023): 71-88.

Nodal discontinuous Galerkin methods for Maxwell's equations in Lorentz-Kerr-Raman medium without nonlinear algebraic solver

The propagation of electromagnetic waves is modeled by time-dependent Maxwell's equations coupled with constitutive laws that describe the responses of the media. In this work, we consider a nonlinear model that describes the electromagnetic wave in an optical medium with the linear Lorentz effect and the cubic nonlinear instantaneous Kerr and delayed Raman effects. Mathematically this model obeys an energy conservative/dissipative law. Though there have been active efforts in designing numerical methods (e.g. of finite difference / finite element /discontinuous Galerkin type) to simulate this model, the methods proposed here are distinctive in that they are free of any nonlinear algebraic solvers. Moreover, in the absence of the Raman effect, our methods also enjoy a provable discrete energy law, and optimal a priori error estimates are further established when the exact solutions are sufficiently smooth. The key ingredients of the new methods include some novel treatment in time discretizations and nodal discontinuous Galerkin spatial discretization for the specific nonlinearities, and they also render a local nature of the methods and hence their suitability for parallel implementation with great efficiency. Numerical experiments are performed to illustrate the accuracy, stability, computational efficiency and parallel scalability of the proposed methods. We further apply the methods to simulate some physically relevant problems in one, two, and three dimensions.

A GPU based accelerated solver for simulation of heat transfer during metal casting process

The metal casting process, which is one of the key drivers of the manufacturing industry, involves several physical phenomena occurring simultaneously like fluid flow, phase change, and heat transfer which affect the casting yield and quality. Casting process modeling involves numerical modeling of these phenomena on a computer. In recent decades, this has become an inevitable tool for foundry engineers to make defect-free castings. To expedite computational time graphics processing units (GPUs) are being increasingly used in the numerical modeling of heat transfer and fluid flow. Initially, in this work a CPU based implicit solver code is developed for solving the 3D unsteady energy equation including phase change numerically using finite volume method which predicts the thermal profile during solidification in the metal casting process in a completely filled mold. To address the computational bottleneck, which is identified as the linear algebraic solver based on the bi-conjugate gradient stabilized method, a GPU-based code is developed using Compute Unified Device Architecture toolkit and was implemented on the GPU. The CPU and GPU based codes are then validated against a commercial casting simulation code FLOW-3D CAST® for a simple casting part and against in-house experimental results for gravity die casting of a simple geometry. Parallel performance is analyzed for grid sizes ranging from 10 × 10 × 10 to 210 × 210 × 210 and for three time-step sizes. The performance of the GPU code based on occupancy and throughput is also investigated. The GPU code exhibits a maximum speedup of 308× compared to the CPU code for a grid size of 210 × 210 × 210 and a time-step size of 2 s.

A multiscale preconditioner for crack evolution in porous microstructures: Accelerating phase‐field methods

AbstractPhase‐field methods are attractive for simulating the mechanical failure of geometrically complex porous microstructures described by 2D/3D x‐ray CT images in subsurface (e.g., CO storage) and manufacturing (e.g., Li‐ion battery) applications. They capture the nucleation, growth, and branching of fractures without prior knowledge of the propagation path or having to remesh the domain. Their drawback lies in the high computational cost for the typical domain sizes encountered in practice. We present a multiscale preconditioner that significantly accelerates the convergence of Krylov solvers in computing solutions of linear(ized) systems arising from the sequential discretization of the momentum and crack‐evolution equations in phase‐field methods. The preconditioner is an algebraic reformulation of a recent pore‐level multiscale method (PLMM) by the authors and consists of a global preconditioner and a local smoother . Together, and attenuate low‐ and high‐frequency errors simultaneously. The proposed , used in the momentum equation only, is a simplification of a recent variant proposed by the authors that is much cheaper and easier to deploy in existing solvers. The smoother , used in both the momentum and crack‐evolution equations, is built such that it is compatible with and more robust and efficient than black‐box smoothers like ILU(). We test and systematically for static‐ and evolving‐crack problems on complex 2D/3D porous microstructures, and show that they outperform existing algebraic multigrid solvers. We also probe different strategies for updating as cracks evolve and show the associated cost can be minimized if is updated adaptively and infrequently. Both and are scalable on parallel machines and can be implemented non‐intrusively in existing codes.

Linear-Time Solution of 3D Magnetostatics With Aggregation-Based Algebraic Multigrid Solvers

Linear-Time Solution of 3D Magnetostatics With Aggregation-Based Algebraic Multigrid Solvers

Light Weight Coarse Grid Aggregation for Smoothed Aggregation Algebraic Multigrid Solver

Light Weight Coarse Grid Aggregation for Smoothed Aggregation Algebraic Multigrid Solver

Evaluation of advection schemes and surface tension model for algebraic and geometric VOF multiphase flow solvers

This paper addresses the validation and verification of numerical models for complex multiphase flows involving immiscible fluids. Specifically, the focus is on the dynamics of the interface between the fluids, which remains a challenge in many applications. Previous studies have primarily examined simpler flows with limited complexity, while real-world scenarios involve intricate interfaces and turbulent dynamics. To overcome this, a representative test case involving the motion of an oil cube within a water cube proposed by [12] is used. Different computational models using volume-of-fluid (VOF) approaches are validated and compared, with particular attention given to the influence of numerical schemes for convective terms. The impact of these schemes on integral quantities and interface dynamics is analyzed. Particular attention is paid at the surface tension contribution to the energy budgets, which is usually neglected in the available literature. Furthermore, it is shown here that, the model adopted for the surface tension is responsible for an additional contribution to numerical dissipation. To the authors' knowledge, it is the first time that a detailed analysis of such aspects is presented.

Unconditionally stable higher order semi-implicit level set method for advection equations

We present compact semi-implicit finite difference schemes on structured grids for numerical solutions of advection by an external velocity and by a speed in the normal direction that are applicable in level set methods. The recommended numerical scheme is third order accurate for the linear advection in the 2D case with a space dependent velocity. Using a combination of analytical and numerical tools in the von Neumann stability analysis, the third order scheme is claimed to be unconditionally stable. We also present a simple high-resolution scheme that gives a TVD (Total Variation Diminishing) approximation of the spatial derivative for the advected level set function in the 1D case. In the case of nonlinear advection, a semi-implicit discretization is proposed to linearize the problem. The compact implicit part of the stencil of numerical schemes contains unknowns only in the upwind direction. Consequently, algebraic solvers like the fast sweeping method can be applied efficiently to solve the resulting algebraic systems. Numerical tests to evolve a smooth and non-smooth interface and an example with a large variation of the velocity confirm the good accuracy of the third order scheme even in the case of very large Courant numbers. The advantage of the high-resolution scheme is documented for examples where the advected level set functions contain large jumps in the gradient.

Agglomeration of polygonal grids using graph neural networks with applications to multigrid solvers

Agglomeration-based strategies are important both within adaptive refinement algorithms and to construct scalable multilevel algebraic solvers. In order to automatically perform agglomeration of polygonal grids, we propose the use of Machine Learning (ML) strategies, that can naturally exploit geometrical information about the mesh in order to preserve the grid quality, enhancing performance of numerical methods and reducing the overall computational cost. In particular, we employ the k-means clustering algorithm and Graph Neural Networks (GNNs) to partition the connectivity graph of a computational mesh. Moreover, GNNs have high online inference speed and the advantage to process naturally and simultaneously both the graph structure of mesh and the geometrical information, such as the areas of the elements or their barycentric coordinates. These techniques are compared with METIS, a standard algorithm for graph partitioning, which is meant to process only the graph information of the mesh. We demonstrate that performance in terms of quality metrics is enhanced for ML strategies. Such models also show a good degree of generalization when applied to more complex geometries, such as brain MRI scans, and the capability of preserving the quality of the grid. The effectiveness of these strategies is demonstrated also when applied to MultiGrid (MG) solvers in a Polygonal Discontinuous Galerkin (PolyDG) framework. In the considered experiments, GNNs show overall the best performance in terms of inference speed, accuracy and flexibility of the approach.

On floating point precision in computational fluid dynamics using OpenFOAM

Thanks to the computational power of modern cluster machines, numerical simulations can provide, with an unprecedented level of details, new insights into fluid mechanics. However, taking full advantage of this hardware remains challenging since data communication remains a significant bottleneck to reaching peak performances. Reducing floating point precision is a simple and effective way to reduce data movement and improve the computational speed of most applications. Nevertheless, special care needs to be taken to ensure the quality and convergence of computed solutions, especially when dealing with complex fluid simulations. In this work, we analyse the impact of reduced (single and mixed compared to double) precision on computational performance and accuracy for computational fluid dynamics. Using the open source library OpenFOAM, we consider incompressible, compressible, and multiphase fluid solvers for testing on relevant benchmarks for flows in the laminar and turbulent regime and in the presence of shock waves. Computational gain and changes in the scalability of applications in reduced precision are also discussed. In particular, an ad hoc theoretical model for the strong scaling allows us to interpret and understand the observed behaviours, as a function of floating point precision and hardware specifics. Finally, we show how reduced precision can significantly speed up a hybrid CPU–GPU implementation, made available to OpenFOAM end-users recently, that simply relies on a GPU linear algebra solver developed by hardware vendors.

Optimization Problems for PDEs in Weak Space-Time Form

Optimization problems constrained by time-dependent Partial Differential Equations (PDEs) are challenging from a computational point of view: even in the simplest case, one needs to solve a system of PDEs coupled globally in time and space for the unknown solutions (the state, the costate and the control of the system). Typical and practically relevant examples are the control of nonlinear heat equations as they appear in laser hardening or the thermic control of flow problems (Boussinesq equations). Specifically for PDEs with a long time horizon, conventional time-stepping methods require an enormous amount of computer memory allocations for the respective other variables. In contrast, adaptive-in-time-and-space methods aim at distributing the available degrees of freedom in an a-posteriori fashion to capture singularities and are, therefore, most promising. Recently, well-posed weak variational formulations have been introduced for time-dependent PDEs such as the heat equation, linear transport and the wave equation. Those formulations also allow for a sharp relation between the approximation error and the residual, which is particularly relevant for model reduction. Moreover, for those tensor-basis formulations, advanced algebraic solvers designed to take into account these multiarray (tensorial) formulations appear to be particularly competitive with respect to time-marching schemes, especially in higher dimensions. We plan to discuss whether these techniques can be extended to nonlinear PDEs like Hamilton-Jacobi-Bellman equations, or stochastic PDEs and variational inequalities. Another topic will be adaptive schemes which, when properly designed, inherit the stability of the continuous formulation. The central goals of the workshop are the analysis, fast solvers and model reduction for PDE-constrained control and optimization problems based on weak formulations of the underlying PDE(s).

Fast large-scale gravity modeling using a high-order compact cascadic multigrid strategy

A new high-order compact (HOC) cascadic multigrid (MG) strategy is developed for modeling large-scale complex gravity problems on nonuniform rectilinear grids. First, the governing 3D Poisson’s gravitational potential equation is discretized by a novel HOC difference scheme, which allows for Robin boundary conditions to further reduce the computational domain. Then, the resulting asymmetric linear system of equations is solved by a generalized extrapolation cascadic MG (gEXCMG) method augmented with a novel MG prolongation operator. The quartic interpolation and completed Richardson extrapolation assist this MG prolongation operator to produce a good initial guess for the biconjugate gradient stabilized smoother. Two synthetic and SEG/EAGE salt models are examined to verify the presented approach. Results indicate that our algorithm can efficiently and accurately produce the gravity signals for complex geologic models. Moreover, the comparisons with the state-of-the-art algebraic MG solvers demonstrate that our gEXCMG solver offers substantially better efficiency with significantly less memory consumption. Therefore, our newly developed method has the potential for serving as a forward engine in large-scale gravity inversions.

Acceleration of supersonic/hypersonic reactive CFD simulations via heterogeneous CPU-GPU supercomputing

The numerical study of reactive flows subjected to supersonic conditions is accelerated by the co-design of a novel strategy to integrate finite-rate chemistry by an adaptive multi-block ODE algebra solver for Graphical Processing Units (GPU), that is coupled to a parallel, shock-capturing Finite-Volume reactive flow solver running on CPUs. The resulting GPGPU solver is validated on Large Eddy Simulations (LES) of a scramjet configuration, whose experimental measurements are available from the literature. It is demonstrated that the proposed method significantly accelerates the solution of reactive CFD computations with Direct Integration of the finite-rate chemistry.

Algebraic Multigrid Using a Stencil–CSR Hybrid Format on GPUs

We propose a new sparse matrix format which captures the matrix structure typical for discretized partial differential equations with piecewise-constant coefficients. The format uses a stencil representation for some blocks of matrix rows, typically corresponding to regions with constant coefficients, whereas other rows are encoded in the general compressed sparse row format. The stencil representation saves memory and is suitable for SIMD-like parallelism as available on GPUs. Further, this format is well suited for the implementation of algebraic multigrid methods, and we present a proof-of-concept GPU-accelerated aggregation-based algebraic multigrid solver based on this format. This solver is compared on a few model problems (2-dimension and 3-dimension Poisson-like) with the compressed sparse row–based solver AmgX from NVIDIA and with the CUDA version of the BoomerAMG solver. For the considered tested problems with one million unknowns or more, the presented solver outperforms AmgX and BoomerAMG in terms of both run time and memory usage, and the performance gap increases with the system size.

Adaptive FEM with quasi-optimal overall cost for nonsymmetric linear elliptic PDEs

Abstract We consider a general nonsymmetric second-order linear elliptic partial differential equation in the framework of the Lax–Milgram lemma. We formulate and analyze an adaptive finite element algorithm with arbitrary polynomial degree that steers the adaptive meshrefinement and the inexact iterative solution of the arising linear systems. More precisely, the iterative solver employs, as an outer loop, the so-called Zarantonello iteration to symmetrize the system and, as an inner loop, a uniformly contractive algebraic solver, for example, an optimally preconditioned conjugate gradient method or an optimal geometric multigrid algorithm. We prove that the proposed inexact adaptive iteratively symmetrized finite element method leads to full linear convergence and, for sufficiently small adaptivity parameters, to optimal convergence rates with respect to the overall computational cost, i.e., the total computational time. Numerical experiments underline the theory.

Controlled-source electromagnetic modeling using a versatile secondary electric-field formulation and efficient multigrid-based preconditioner

Frequency-domain controlled-source electromagnetic (CSEM) surveying is a widely used geophysical electromagnetic exploration method. To avoid the singularity near the transmitting sources and improve the accuracy of three-dimensional (3D) CSEM forward modeling, we adopt the vector finite element method and secondary electric field formulation combined with hexahedral grids to solve the frequency-domain Maxwell equations. Primary field computation applies the fast Hankel transform using digital linear filters to quickly calculate the electromagnetic field generated by different electromagnetic dipoles or loops in the layered medium. Therefore, our forward modeling algorithm is highly versatile and compatible with various electromagnetic sources and scenarios of CSEM application (e.g., airborne, offshore, or marine). To improve the convergence rate, we use an algebraic multigrid solver based on the established Hiptmair–Xu decomposition. Using direct solvers requires a vast amount of memory, while conventional iterative solvers usually fail or converge only very slowly. Our new iterative scheme does not require the divergence corrections commonly applied in finite difference or finite volume electromagnetic modeling to accelerate the convergence. We implement complete parallelism in our finite element programming to efficiently solve the CSEM forward problems on massively parallel clusters. We use different numerical methods to obtain high-precision comparison results for one-dimensional and typical 3D models to demonstrate our algorithm's effectiveness and correctness. Finally, a public model of 3 million cells is used to simulate the responses of marine CSEM surveying. The proposed new solver shows robust scalability when using nearly 1000 cores on the cluster, demonstrating that our algorithm is particularly suitable for large-scale 3D forward modeling.

A Dimension-Oblivious Domain Decomposition Method Based on Space-Filling Curves

In this paper we present an algebraic dimension-oblivious two-level domain decomposition solver for discretizations of elliptic partial differential equations. The proposed parallel domain decomposition solver is based on a space-filling curve partitioning approach that is applicable to any discretization, i.e., it directly operates on the assembled matrix equations. Moreover, it allows for the effective use of arbitrary processor numbers independent of the dimension of the underlying partial differential equation while maintaining its optimal convergence behavior. This is the core property required to attain a combination technique solver with extreme scalability which can utilize exascale parallel systems efficiently. Moreover, this approach provides a basis for the development of a fault-tolerant solver for the numerical treatment of high-dimensional Problems. To achieve the required data redundancy we are therefore concerned with large overlaps of our domain decomposition, which we construct in any dimension via space-filling curves. In this paper, we propose a space-filling curve based domain decompositoin solver and present its convergence properties and scaling behavior. The results of our numerical experiments clearly show that our approach provides optimal convergence and scaling behavior in an arbitrary dimension utilizing arbitrary processor numbers.

Guaranteed contraction of adaptive inexact hp-refinement strategies with realistic stopping criteria

The purpose of this contribution it to theoretically analyze the adaptive refinement strategies for conforming hp-finite element approximations of elliptic problems proposed for exact algebraic solvers in (P. Daniel, A. Ern, I. Smears and M. Vohralk, Comput. Math. Appl. 76 (2018) 967–983.) and for inexact algebraic solvers in (P. Daniel, A. Ern, and M. Vohralk, Comput. Methods Appl. Mech. Eng. 359 (2020) 112607.). Both of these strategies are driven by guaranteed equilibrated flux energy error estimators. The employed hp-refinement criterion stems from solving two separate local residual problems posed only on the patches of elements around marked vertices selected by a bulk-chasing criterion. In the above references, we have derived a fully computable guaranteed bound on the ratio of the error on two successive steps of the hp-adaptive loop. Here, our focus is to prove that this ratio is uniformly smaller than one, and thus the convergence of the adaptive and adaptive inexact hp-refinement strategies. To be able to achieve this goal, we have to introduce some additional assumptions on the h- and p-refinements, namely an extension of the marked region, as well as a sufficient h- or p-refinement of each marked patch. We investigate two such strategies, where one ensures a polynomial-degree-robust guaranteed contraction. In the inexact case, a sufficiently precise stopping criterion for the algebraic solver is requested, but this criterion remains fully computable and also realistic in the sense that in our numerical experiments, it does not request the algebraic error to be excessively small in comparison with the total error.

Stationary Stokes solver for single-phase flow in porous media: A blastingly fast solution based on Algebraic Multigrid Method using GPU

The paper is focused on high efficiency Stokes solver that is applied to the incompressible flow in porous media. Computational domains are represented by binarized 3D computed tomography voxel models. The solution procedure is constructed around a fast algebraic multigrid solver that utilizes the power of graphics processing units (GPUs). In order to minimize memory footprint and accelerate the solver, a simple MAC-type staggered finite difference discretization is used and a coupled Stokes saddle-point type system is solved directly on a GPU. The MAC discretization is discussed, taking particular care of internal solid boundaries around pores. The method includes topological domain analysis on a GPU, which removes isolated volumes with no flow in the domain (based on the connected component labeling), depending on the boundary conditions. We consider various types of boundary conditions and efficient parallel strategies for the GPUs, including fast matrix assembly and residual regularization. Our method is extensively benchmarked both against analytical solutions and applied problems from digital rock physics in terms of computational wall time, precision, approximation order, and convergence. We demonstrate that it takes up to 5–23 s on a modern GPU card to obtain a solution with 1⋅10−6 error residuals for 3D geometries with 300–4503 voxels and a porosity range of 5–37%.

An efficient extrapolation multigrid method based on a HOC scheme on nonuniform rectilinear grids for solving 3D anisotropic convection–diffusion problems

We develop an efficient multigrid method combined with a high-order compact (HOC) finite difference scheme on nonuniform rectilinear grids for solving 3D diagonal anisotropic convection–diffusion problems with boundary/interior layers. Firstly, we derive a fourth-order compact finite difference scheme to discretize the 3D anisotropic convection–diffusion equation on a rectilinear grid. Then, the resulting large-scale asymmetric linear system of equations is solved by a generalized extrapolation cascadic multigrid (gEXCMG) method based on two novel multigrid (MG) prolongation operators. The highlight of this paper is the application of the quintic Lagrange interpolation and the completed Richardson extrapolation in the design of the new MG prolongation operator on nonuniform rectilinear grids, which can produce a good initial guess (sixth-order approximation to the finite difference solution) for the SSOR-preconditioned biconjugate gradient stabilized (BiCGStab) smoother. In the end, numerical experiments show that the gEXCMG method combined with the HOC scheme can achieve fourth-order accuracy for 3D anisotropic convection–diffusion problems with few smoothing steps on the finest grid. Moreover, the proposed gEXCMG method can offer substantially better efficiency than the state-of-the-art algebraic MG solver, namely, aggregation-based algebraic multigrid (AGMG) method, for large linear systems arising from the discretization of second order elliptic PDEs.