Stand-alone Solver Research Articles

AI‐enhanced iterative solvers for accelerating the solution of large‐scale parametrized systems

SummaryRecent advances in the field of machine learning open a new era in high performance computing for challenging computational science and engineering applications. In this framework, the use of advanced machine learning algorithms for the development of accurate and cost‐efficient surrogate models of complex physical processes has already attracted major attention from scientists. However, despite their powerful approximation capabilities, surrogate model predictions are still far from being near to the ‘exact’ solution of the problem. To address this issue, the present work proposes the use of up‐to‐date machine learning tools in order to equip a new generation of iterative solvers of linear equation systems, capable of very efficiently solving large‐scale parametrized problems at any desired level of accuracy. The proposed approach consists of the following two steps. At first, a reduced set of model evaluations is performed using a standard finite element methodology and the corresponding solutions are used to establish an approximate mapping from the problem's parametric space to its solution space using a combination of deep feedforward neural networks and convolutional autoencoders. This mapping serves as a means of obtaining very accurate initial predictions of the system's response to new query points at negligible computational cost. Subsequently, an iterative solver inspired by the Algebraic Multigrid method in combination with Proper Orthogonal Decomposition, termed POD‐2G, is developed that successively refines the initial predictions of the surrogate model towards the exact solution. The application of POD‐2G as a standalone solver or as preconditioner in the context of preconditioned conjugate gradient methods is demonstrated on several numerical examples of large scale systems, with the results indicating its strong superiority over conventional iterative solution schemes.

MGM: A Meshfree Geometric Multilevel Method for Systems Arising from Elliptic Equations on Point Cloud Surfaces

We develop a new meshfree geometric multilevel (MGM) method for solving linear systems that arise from discretizing elliptic PDEs on surfaces represented by point clouds. The method uses a Poisson disk sampling-type technique for coarsening the point clouds and new meshfree restriction/interpolation operators based on polyharmonic splines for transferring information between the coarsened point clouds. These are then combined with standard smoothing and operator coarsening methods in a V-cycle iteration. MGM is applicable to discretizations of elliptic PDEs based on various localized meshfree methods, including RBF finite differences (RBF-FD) and generalized finite differences (GFD). We test MGM both as a standalone solver and preconditioner for Krylov subspace methods on several test problems using RBF-FD and GFD and numerically analyze convergence rates, efficiency, and scaling with increasing point cloud sizes. We also perform a side-by-side comparison to algebraic multigrid methods for solving the same systems. Finally, we further demonstrate the effectiveness of MGM by applying it to three challenging applications on complicated surfaces: pattern formation, surface harmonics, and geodesic distance.

On the Recursive Structure of Multigrid Cycles

A new non-adaptive recursive scheme for multigrid algorithms is introduced. Governed by a positive parameter $\kappa$ called the cycle counter, this scheme generates a family of multigrid cycles dubbed $\kappa$-cycles. The well-known $V$-cycle, $F$-cycle, and $W$-cycle are shown to be particular members of this rich $\kappa$-cycle family, which satisfies the property that the total number of recursive calls in a single cycle is a polynomial of degree $\kappa$ in the number of levels of the cycle. This broadening of the scope of fixed multigrid cycles is shown to be potentially significant for the solution of some large problems on platforms, such as GPU processors, where the overhead induced by recursive calls may be relatively significant. In cases of problems for which the convergence of standard $V$-cycles or $F$-cycles (corresponding to $\kappa=1$ and $\kappa=2$, respectively) is particularly slow, and yet the cost of $W$-cycles is very high due to the large number of recursive calls (which is in exponential in the number of levels), intermediate values of $\kappa$ may prove to yield significantly faster run-times. This is demonstrated in examples where $\kappa$-cycles are used for the solution of rotated anisotropic diffusion problems, both as a stand-alone solver and as a preconditioner. Moreover, a simple model is presented for predicting the approximate run-time of the $\kappa$-cycle, which is useful in pre-selecting an appropriate cycle counter for a given problem on a given platform. Implementing the $\kappa$-cycle requires making just a small change in the classical multigrid cycle.

Machine Learning Assisted Stochastic Unit Commitment During Hurricanes With Predictable Line Outages

Stochastic unit commitment is an efficient method for grid operation in the presence of significant uncertainties. An example is an operation during a predicted hurricane with uncertain line out-ages. However, the solution quality comes at the cost of substantial computational burden, which makes its adoption challenging. This paper evaluates some possible ways that machine learning can be used to reduce this computational burden. First, a set of feasibility studies is conducted. Results suggest that using machine learning as an assistant to the stochastic unit commitment solver is more advantageous than using it as a standalone solver. In particular, the machine learning model is trained to facilitate solving the problem by determining the unnecessary constraints that can be removed from the original problem without affecting the final accuracy. The variables that can be used as input features/predictors or outputs for the machine learning model are determined through feasibility studies. Then, an algorithm to train and utilize a machine learning model is proposed. The method is tested on a 500-bus synthetic South Carolina system. Various test cases show an average reduction in solution time by more than 90% by using the trained machine learning model to assist the stochastic unit commitment solver.

ASUSD nuclear data sensitivity and uncertainty program package: Validation on fusion and fission benchmark experiments

Nuclear data (ND) sensitivity and uncertainty (S/U) quantification in shielding applications is performed using deterministic and probabilistic approaches. In this paper the validation of the newly developed deterministic program package ASUSD (ADVANTG + SUSD3D) is presented. ASUSD was developed with the aim of automating the process of ND S/U while retaining the computational efficiency of the deterministic approach to ND S/U analysis. The paper includes a detailed description of each of the programs contained within ASUSD, the computational workflow and validation results.ASUSD was validated on two shielding benchmark experiments from the Shielding Integral Benchmark Archive and Database (SINBAD) - the fission relevant ASPIS Iron 88 experiment and the fusion relevant Frascati Neutron Generator (FNG) Helium Cooled Pebble Bed (HCPB) Test Blanket Module (TBM) mock-up experiment. The validation process was performed in two stages. Firstly, the Denovo discrete ordinates transport solver was validated as a standalone solver. Secondly, the ASUSD program package as a whole was validated as a ND S/U analysis tool. Both stages of the validation process yielded excellent results, with a maximum difference of 17% in final uncertainties due to ND between ASUSD and the stochastic ND S/U approach. Based on these results, ASUSD has proven to be a user friendly and computationally efficient tool for deterministic ND S/U analysis of shielding geometries.

SCALAR: an AMR code to simulate axion-like dark matter models

We present a new code, SCALAR, based on the high-resolution hydrodynamics and N-body code RAMSES, to solve the Schrödinger equation on adaptive refined meshes. The code is intended to be used to simulate axion or fuzzy dark matter models where the evolution of the dark matter component is determined by a coupled Schrödinger-Poisson equation, but it can also be used as a stand-alone solver for both linear and non-linear Schrödinger equations with any given external potential. This paper describes the numerical implementation of our solver and presents tests to demonstrate how accurately it operates.

Solver Requirements for Interactive Configuration

Interactive configuration includes the user as an essential factor in the configuration process. The two main components of an interactive configurator are a user interface at the front-end and a knowledge representation and reasoning (KRR) framework at the back-end. In this paper we discuss important requirements for the underlying KRR system to support an interactive configuration process. Representative of many reasoning systems and tools used for implementing product configurators, we selected MiniZinc, Choco, Potassco, Picat, CP-SAT solver, and Z3 for evaluation and reviewed them against the identified requirements. We observe that many of those requirements are not well supported by existing stand-alone solvers.

A computational investigation of preconditioning strategies and iterative methods for finite element based neurostimulation simulations

Computational simulations of transcranial direct current stimulation (tDCS) enable researchers and medical practitioners to investigate this form of neurostimulation with in silico experiments. For these computer-based simulations to be of practical use to the medical community, patient-specific head geometries and finely discretized computational grids must be used. As a result, solving the partial differential equation based mathematical model that governs tDCS can be computationally burdensome. Further, consider the task of identifying optimal electrode configurations and parameters for a particular patient’s condition, head geometry, and therapeutic objectives; it is conceivable that hundreds of tDCS simulations could be executed. It is therefore important and necessary to identify efficient solution methods for medically-based tDCS simulations. To address this requirement, we exhaustively compare the convergence performance of geometric multigrid and the preconditioned conjugate gradient method when solving the linear system of equations generated from a finite element discretization of the tDCS governing equations. Our simulations consist of three commonly used real-world tDCS electrode montages on MRI-derived three-dimensional head models with physiologically-based inhomogeneous tissue conductivities. Simulations are realized on fine computational meshes with resolutions deemed applicable to the medical community, and as a result, our finite element solutions highlight tDCS-specific phenomena such as electric field shunting that contributes to a notable intensification of the stimulation’s electric current dosage. Convergence metrics of each linear system solver are examined, and compared and linked to theoretical estimates. It is shown that the conjugate gradient method achieves superior convergence rates only when preconditioned with an appropriately configured multigrid algorithm. In addition, it is demonstrated that physiological characteristics of tDCS simulations make multigrid as a stand-alone solver highly ineffective, despite the fact that this method is typically effective in solving the tDCS-based Poisson problem. By identifying the solution methods optimal for medically-driven tDCS simulations, our results extend simulation support to high-resolution and high-volume computing applications, and will ultimately help guide tDCS numerical simulations towards becoming an integrated aspect of the patient-specific tDCS treatment protocol.

Outer approximation for integer nonlinear programs via decision diagrams

As an alternative to traditional integer programming (IP), decision diagrams (DDs) provide a new solution technology for discrete problems exploiting their combinatorial structure and dynamic programming representation. While the literature mainly focuses on the competitive aspects of DDs as a stand-alone solver, we investigate their complementary role by introducing IP techniques that can be derived from DDs and used in conjunction with IP to enhance the overall performance. This perspective allows for studying problems with more general structure than those typically modeled via recursive formulations. In particular, we develop linear programming and subgradient-type methods to generate valid inequalities for the convex hull of the feasible region described by DDs. For convex IPs, these cutting planes dominate the so-called linearized cuts used in the outer approximation framework. These cutting planes can also be derived for nonconvex IPs, which leads to a generalization of the outer approximation framework. Computational experiments show significant optimality gap improvement for integer nonlinear programs over the traditional cutting plane methods employed in the state-of-the-art solvers.

Modeling and Optimization of Truck-Shovel Allocation to Mining Faces in Cement Quarry

Truck and shovel are the most common raw material transportation system used in the cement quarry operations. One of the major challenges associated with the cement quarry operations is the efficient allocation of truck and shovel to the mining faces. In order to minimize the truck and shovel operating cost, subject to quantity and quality constraints, the mixed integer linear programing (MILP) model for truck and shovel allocation to mining faces for cement quarry is presented. This model is implemented using the optimization IDE tool GUSEK (GLPK under SciTE Extended Kit) and the GLPK (GNU Linear Programming Kit) standalone solver. The MILP model is applied to an existing cement quarry operation, the Kohat cement quarry located at Kohat (Pakistan) as a case study. The analysis of the results of the relating case study reveals that significant gains are achievable through employing the MILP model. The results obtained not only show a significant cost reduction but also help in achieving a better coordination among the quarry and quality department.

Scheduled relaxation Jacobi method as preconditioner of Krylov subspace techniques for large-scale Poisson problems

This article presents an assessment of the scheduled relaxation Jacobi (SRJ) method for the solution of large-scale Poisson problems arising in the numerical simulation of large eddy turbulent flow in large complex geometry. The SRJ schemes are used both as standalone solvers and as preconditioners to Krylov subspace solvers. The Navier-Stokes equation is solved using structured Cartesian grids (both uniform as well as nonuniform grids). Geometrical complexities are handled using the immersed boundary method (IBM) method. The performance of SRJ schemes as a standalone solver is first validated with the help of a rectangular channel flow problem whose exact solution is known. Numerical experiments are performed to evaluate the performance of SRJ schemes as preconditioners to two robust and efficient Krylov subspace solvers (PCG and BiCGSTAB), and compared with different preconditioners such as Jacobi (diagonal preconditioner), SOR(k) and multigrid preconditioners. SRJ schemes as standalone solver perform quite similarly as predicted before in literature. However, as preconditioner particularly some SRJ schemes work best among the rest for three-dimensional problems.

A Distributed Approach to LARS Stream Reasoning (System paper)

AbstractStream reasoning systems are designed for complex decision-making from possibly infinite, dynamic streams of data. Modern approaches to stream reasoning are usually performing their computations using stand-alone solvers, which incrementally update their internal state and return results as the new portions of data streams are pushed. However, the performance of such approaches degrades quickly as the rates of the input data and the complexity of decision problems are growing. This problem was already recognized in the area of stream processing, where systems became distributed in order to allocate vast computing resources provided by clouds. In this paper we propose a distributed approach to stream reasoning that can efficiently split computations among different solvers communicating their results over data streams. Moreover, in order to increase the throughput of the distributed system, we suggest an interval-based semantics for the LARS language, which enables significant reductions of network traffic. Performed evaluations indicate that the distributed stream reasoning significantly outperforms existing stand-alone LARS solvers when the complexity of decision problems and the rate of incoming data are increasing.

On a preconditioner for time domain boundary element methods

We propose a time stepping scheme for the space-time systems obtained from Galerkin time-domain boundary element methods for the wave equation. Based on extrapolation, the method proves stable, becomes exact for increasing degrees of freedom and can be used either as a preconditioner, or as an efficient standalone solver for scattering problems with smooth solutions. It also significantly reduces the number of GMRES iterations for screen problems, with less regularity, and we explore its limitations for enriched methods based on non-polynomial approximation spaces.

BootCMatch

This article has two main objectives: one is to describe some extensions of an adaptive Algebraic Multigrid (AMG) method of the form previously proposed by the first and third authors, and a second one is to present a new software framework, named BootCMatch , which implements all the components needed to build and apply the described adaptive AMG both as a stand-alone solver and as a preconditioner in a Krylov method. The adaptive AMG presented is meant to handle general symmetric and positive definite (SPD) sparse linear systems, without assuming any a priori information of the problem and its origin; the goal of adaptivity is to achieve a method with a prescribed convergence rate. The presented method exploits a general coarsening process based on aggregation of unknowns, obtained by a maximum weight matching in the adjacency graph of the system matrix. More specifically, a maximum product matching is employed to define an effective smoother subspace (complementary to the coarse space), a process referred to as compatible relaxation, at every level of the recursive two-level hierarchical AMG process. Results on a large variety of test cases and comparisons with related work demonstrate the reliability and efficiency of the method and of the software.

On preconditioned BiCGSTAB solver for MLPG method applied to heat conduction in 3D complex geometry

Development of meshfree methods is in full swing in the last few decades as these methods overcome the shortcomings of mesh based methods. Development of meshfree methods for large scale problems requires fast iterative solvers. Preconditioner plays an important role in improving convergence rate of an iterative solver. Recently proposed SRJ method can be used as a stand-alone solver and as a preconditioner in CG type methods. Objective of the present work is to evaluate the performance of SRJ method as a preconditioner in BiCGSTAB method in connection with the MLPG method. Three dimensional Poisson equation and heat conduction equation are solved in regular and irregular domain respectively. The SRJ method is tested as a solver and as a preconditioner in BiCGSTAB method. Its performance as a preconditioner is compared with Jacobi and SOR preconditioners.

On preconditioned BiCGSTAB solver for MLPG method applied to heat conduction in complex geometry

ABSTRACTMeshless local Petrov–Galerkin (MLPG) method is a promising meshfree method for continuum problems in complex domains, especially for large deformation, moving boundary and phase change problems. For large-scale problems, iterative methods for solving the discretized equations are more suitable than direct methods. Krylov subspace solvers of conjugate gradient type are the most preferred iterative solvers. The convergence rate of these methods depends on preconditioner used. Recently, proposed schedule relaxation Jacobi (SRJ) method can be used as a stand-alone solver and as a preconditioner. In the present work, the SRJ method is tested as a stand-alone solver and as a preconditioner for BiCGSTAB solver in the MLPG method, and its performance has been compared with successive overrelaxation (k) preconditioner. Two-dimensional linear steady-state heat conduction in complex shape geometry has been used as the model test problem.

Hybrid Fourier pseudospectral/discontinuous Galerkin time-domain method for wave propagation

The Fourier Pseudospectral time-domain (Fourier PSTD) method was shown to be an efficient way of modelling acoustic propagation problems as described by the linearized Euler equations (LEE), but is limited to real-valued frequency independent boundary conditions and predominantly staircase-like boundary shapes. This paper presents a hybrid approach to solve the LEE, coupling Fourier PSTD with a nodal Discontinuous Galerkin (DG) method. DG exhibits almost no restrictions with respect to geometrical complexity or boundary conditions. The aim of this novel method is to allow the computation of complex geometries and to be a step towards the implementation of frequency dependent boundary conditions by using the benefits of DG at the boundaries, while keeping the efficient Fourier PSTD in the bulk of the domain. The hybridization approach is based on conformal meshes to avoid spatial interpolation of the DG solutions when transferring values from DG to Fourier PSTD, while the data transfer from Fourier PSTD to DG is done utilizing spectral interpolation of the Fourier PSTD solutions. The accuracy of the hybrid approach is presented for one- and two-dimensional acoustic problems and the main sources of error are investigated. It is concluded that the hybrid methodology does not introduce significant errors compared to the Fourier PSTD stand-alone solver. An example of a cylinder scattering problem is presented and accurate results have been obtained when using the proposed approach. Finally, no instabilities were found during long-time calculation using the current hybrid methodology on a two-dimensional domain.

Computational fluid dynamics for nematic liquid crystals

Due to recent advances in fast iterative solvers in the field of computational fluid dynamics, more complex problems which were previously beyond the scope of standard techniques can be tackled. In this paper, we describe one such situation, namely, modelling the interaction of flow and molecular orientation in a complex fluid such as a liquid crystal. Specifically, we consider a nematic liquid crystal in a spatially inhomogeneous flow situation where the orientational order is described by a second rank alignment tensor. The evolution is determined by two coupled equations: a generalised Navier–Stokes equation for flow in which the divergence of the stress tensor also depends on the alignment tensor and its time derivative, and a convection-diffusion type equation with non-linear terms that stem from a Landau-Ginzburg-DeGennes potential for the alignment. In this paper, we use a specific model with three viscosity coefficients that allows the contribution of the orientation to the viscous stress to be cast in the form of an orientation-dependent force. This effectively decouples the flow and orientation, with each appearing only on the right-hand side of the other equation. In this way, difficulties associated with solving the fully coupled problem are circumvented and a stand-alone fast solver, such as the state-of-the-art preconditioned iterative solver implemented here, can be used for the flow equation. A time-discretised strategy for solving the flow-orientation problem is illustrated using the example of Stokes flow in a lid-driven cavity.

A multiscale two-point flux-approximation method

A large number of multiscale finite-volume methods have been developed over the past decade to compute conservative approximations to multiphase flow problems in heterogeneous porous media. In particular, several iterative and algebraic multiscale frameworks that seek to reduce the fine-scale residual towards machine precision have been presented. Common for all such methods is that they rely on a compatible primal–dual coarse partition, which makes it challenging to extend them to stratigraphic and unstructured grids. Herein, we propose a general idea for how one can formulate multiscale finite-volume methods using only a primal coarse partition. To this end, we use two key ingredients that are computed numerically: (i) elementary functions that correspond to flow solutions used in transmissibility upscaling, and (ii) partition-of-unity functions used to combine elementary functions into basis functions. We exemplify the idea by deriving a multiscale two-point flux-approximation (MsTPFA) method, which is robust with regards to strong heterogeneities in the permeability field and can easily handle general grids with unstructured fine- and coarse-scale connections. The method can easily be adapted to arbitrary levels of coarsening, and can be used both as a standalone solver and as a preconditioner. Several numerical experiments are presented to demonstrate that the MsTPFA method can be used to solve elliptic pressure problems on a wide variety of geological models in a robust and efficient manner.

Groundwater modelling: Towards an estimation of the acceleration factors of iterative methods via an analysis of the transmissivity spatial variability

When running a groundwater flow model, a recurrent and seemingly subsidiary question arises at the starting step of computations: what value of acceleration parameter do we need to optimize the numerical solver? A method is proposed to provide a practical estimate of the optimal acceleration parameter via a geostatistical analysis of the spatial variability of the logarithm of the transmissivity field Y. The background of the approach is illustrated on the successive over-relaxation method (SOR) used, either as a stand-alone solver, or as a symmetric preconditioner (SSOR) to the gradient conjugate method, or as a smoother in multigrid methods. It shows that this optimum acceleration factor is a function of the standard deviation and the correlation length of Y. This provides an easy-to-use heuristic procedure to estimate the acceleration factors, which could even be incorporated in the software package. A case study illustrates the steps needed to perform this estimation.