Robustness Of Solvers Research Articles

An entropy conserving/stable discontinuous Galerkin solver in entropy variables based on the direct enforcement of entropy balance

Numerical schemes that guarantee the preservation of flow properties as prescribed by physical laws, ensure high-fidelity results and improved computational robustness. This work presents a highly accurate entropy conserving/stable solver for the solution of the compressible Euler equations. The method uses a modal Discontinuous Galerkin approximation and takes advantage of both entropy variables and the Direct Enforcement of Entropy Balance (DEEB), proposed by Abgrall [1] to fulfill the second law of thermodynamics at the discrete level. Thanks to the use of an orthonormal modal basis, the cost of adding the DEEB correction to the spatially discretized equations is negligible. Entropy variables are used directly, as unknowns, or to assemble the discrete operators while solving for the conservative variables. In the latter case, the conservative variables are “converted” into the entropy ones through an L2 “entropy projection”. If this strategy is coupled with DEEB, a very efficient alternative to the direct use of the entropy variables is obtained, with a dramatic improvement in robustness over the direct use of the conservative variables. For the time integration both an explicit Strong Stability Preserving Runge-Kutta scheme and the implicit Generalized Crank-Nicolson method are considered. The latter allows for the construction of a discrete entropy conservative/stable scheme both in space and time to be used as a reference. Convergence and conservation properties are demonstrated together with robustness and computational efficiency on a suite of two-dimensional test cases, showing that by coupling DEEB with entropy projection the solver robustness is improved while not affecting the stability limit.

Towards Fast and Robust Simulation in Integrated Electricity and Gas System: A Sequential United Method

The decoupling method, in which different energy vectors are solved alternatively, is widely adopted in the integrated electricity and gas systems (IEGS) analysis. Nevertheless, the bidirectional interaction created by power to gas (P2G) and gas turbines (GT) requires iterative calculations between the electric power system (EPS) and the natural gas network (NGN). Without proper solutions, the iterations will raise predictable concerns about simulation efficiency and convergence. Therefore, this paper investigates the fast and robust simulation in IEGS considering the nonlinearity and dynamics from both the network and device sides. To address the bidirectional interaction, a state decoupling method (SDM) is first proposed as an efficient solver for the dynamic NGN model. Based on this, a sequential united method (SUM) is developed to reduce the computational cost and improve the solver robustness, which decouples the large-scale Jacobian matrix inversion into the sequential solution of multiple small sub-matrixes. Case studies in two systems demonstrate the superiority of the proposed method compared with the decoupling method, indicating a significant improvement of SUM in simulation efficiency and robustness.

Efficient solution of bimaterial Riemann problems for compressible multi-material flow simulations

When solving compressible multi-material flow problems, an unresolved challenge is the computation of advective fluxes across material interfaces that separate drastically different thermodynamic states and relations. A popular idea in this regard is to locally construct bimaterial Riemann problems, and to apply their exact solutions in flux computation. For general equations of state, however, finding the exact solution of a Riemann problem is expensive as it requires nested loops. Multiplied by the large number of Riemann problems constructed during a simulation, the computational cost often becomes prohibitive. The work presented in this paper aims to accelerate the solution of bimaterial Riemann problems without introducing approximations or offline precomputation tasks. The basic idea is to exploit some special properties of the Riemann problem equations, and to recycle previous solutions as much as possible. Following this idea, four acceleration methods are developed, including (1) a change of integration variable through rarefaction fans, (2) storing and reusing integration trajectory data, (3) step size adaptation, and (4) constructing an R-tree on the fly to generate initial guesses. The performance of these acceleration methods is assessed using four example problems in underwater explosion, laser-induced cavitation, and hypervelocity impact. These problems exhibit strong shock waves, large interface deformation, contact of multiple (>2) interfaces, and interaction between gases and condensed matters. For all the problems, the acceleration methods are able to significantly reduce the computational cost without affecting solver robustness or solution accuracy. In different cases, the solution of bimaterial Riemann problems is accelerated by 37 to 87 times. As a result, the total cost of advective flux computation, which includes the exact Riemann problem solution at material interfaces and the numerical flux calculation over the entire computational domain, is accelerated by 18 to 81 times.

An effective preconditioning strategy for volume penalized incompressible/low Mach multiphase flow solvers

The volume penalization (VP) or the Brinkman penalization (BP) method is a diffuse interface method for simulating multiphase fluid-structure interaction (FSI) problems in ocean engineering and/or phase change problems in thermal sciences and engineering. The method relies on a penalty factor (which is inversely related to body's permeability κ) that must be large to enforce rigid body velocity in the solid domain. When the penalty factor is large, the discrete system of equations becomes stiff and difficult to solve numerically. In this paper, we propose a projection method-based preconditioning strategy for solving volume penalized (VP) incompressible and low-Mach Navier-Stokes equations. The projection preconditioner enables the monolithic solution of the coupled velocity-pressure system in both single phase (uniform density and viscosity) and multiphase (variable density and viscosity) flow settings. In this approach, the penalty force is treated implicitly, which is allowed to take arbitrary large values without affecting the solver's convergence rate or causing numerical stiffness. It is made possible by including the penalty term in the pressure Poisson equation (PPE), which was not included in previous works that solved VP incompressible Navier-Stokes equations using the projection method. We show how and where the Brinkman penalty term enters the PPE by re-deriving the projection algorithm for the VP method. Solver scalability under grid refinement is demonstrated, i.e., convergence is achieved with the same number of iterations regardless of the problem size. A manufactured solution in a single phase setting is used to determine the spatial accuracy of the penalized solution. Various values of body's permeability κ are considered. Second-order pointwise accuracy is achieved for both velocity and pressure solutions for reasonably small values of κ. Error saturation occurs when κ is extremely small, but the convergence rate of the solver does not degrade. The solver converges faster as κ decreases, contrary to prior experience. Two multiphase fluid-structure interaction (FSI) problems from the ocean engineering literature are also simulated to evaluate the solver's robustness and performance (in terms of its number of iterations). The proposed solver also allows us to investigate the effect of κ on the motion of the contact line over the surface of the immersed body. It also allows us to investigate the dynamics of the free surface of a solidifying metal.

Scalable Graphics Processing Unit–Based Multiscale Linear Solvers for Reservoir Simulation

Summary In this work, the scalability of two key multiscale solvers for the pressure equation arising from incompressible flow in heterogeneous porous media, namely, the multiscale finite volume (MSFV) solver, and the restriction-smoothed basis multiscale (MsRSB) solver, are investigated on the graphics processing unit (GPU) massively parallel architecture. The robustness and scalability of both solvers are compared against their corresponding carefully optimized implementation on the shared-memory multicore architecture in a structured problem setting. Although several components in MSFV and MsRSB algorithms are directly parallelizable, their scalability on the GPU architecture depends heavily on the underlying algorithmic details and data-structure design of every step, where one needs to ensure favorable control and data flow on the GPU, while extracting enough parallel work for a massively parallel environment. In addition, the type of algorithm chosen for each step greatly influences the overall robustness of the solver. Thus, we extend the work on the parallel multiscale methods of Manea et al. (2016) to map the MSFV and MsRSB special kernels to the massively parallel GPU architecture. The scalability of our optimized parallel MSFV and MsRSB GPU implementations are demonstrated using highly heterogeneous structured 3D problems derived from the SPE10 Benchmark (Christie and Blunt 2001). Those problems range in size from millions to tens of millions of cells. For both solvers, the multicore implementations are benchmarked on a shared-memory multicore architecture consisting of two packages of Intel® Cascade Lake Xeon Gold 6246 central processing unit (CPU), whereas the GPU implementations are benchmarked on a massively parallel architecture consisting of NVIDIA Volta V100 GPUs. We compare the multicore implementations to the GPU implementations for both the setup and solution stages. Finally, we compare the parallel MsRSB scalability to the scalability of MSFV on the multicore (Manea et al. 2016) and GPU architectures. To the best of our knowledge, this is the first parallel implementation and demonstration of these versatile multiscale solvers on the GPU architecture. NOTE: This paper is also published as part of the 2021 SPE Reservoir Simulation Conference Special Issue.

Computational fluid-dynamic modelling of two-phase compressible flows of carbon dioxide in supercritical conditions

Abstract Compressible two-phase flows of carbon dioxide in supercritical thermodynamic conditions are encountered in many applications, e.g. ejectors for refrigeration and compressors for power production and carbon capture and sequestration to name a few. Alongside the phase change, transonic/supersonic flow regimes and non-ideal effects also add additional complexities in the simulations of such flows. In this work, we investigate cavitating and condensing flows of carbon dioxide via numerical simulations based on the two-fluid concept, applying both a mixture model and a barotropic model. In the mixture model, the phase change is modelled with an extra transport equation for the mass of the dispersed phase and a source term introduced via a penalty formulation. The barotropic model reproduces the pressure–density relation of the mixture along the upstream isentrope. Both the models assume thermodynamic and mechanical equilibrium between phases and exclude meta-stability effects. All results are compared against experimental data taken from literature and the main numerical issues of the models are discussed in detail. The agreement between the simulations and the experiments is remarkable qualitatively and quantitatively, resulting in the range 2%–4% for pressure and below 1% for temperature in terms of weighted mean absolute percentage error for supercritical expansions, even though suggesting a further margin of improvement in the physical modelling, especially for subcritical expansions. Finally, we show that the barotropic model yields comparable predictions of the expansion processes at a lower computational cost and with an improved solver robustness.

ADflow: An Open-Source Computational Fluid Dynamics Solver for Aerodynamic and Multidisciplinary Optimization

Computational fluid dynamics through the solution of the Navier–Stokes equations with turbulence models has become commonplace. However, simply solving these equations is not sufficient to be able to perform efficient design optimization with a flow solver in the loop. This paper discusses the recommendations for developing a flow solver that is suitable for efficient aerodynamic and multidisciplinary design optimization. One of the major recommendations is to be able to load the flow solver as a library that provides direct memory access to the relevant data. Other recommendations are to use a higher-level language for scripting and to pay special attention to solution warm starting, code efficiency, flow solver robustness, and solution failure handling. As an example of a flow solver that follows these recommendations, the open-source flow solver ADflow is presented. Results from aerodynamic optimization, aerostructural analysis, and aerostructural optimization using ADflow demonstrate the performance advantages claimed in the recommendations. The publication of these recommendations and the availability of the source code open the door for other solvers to adopt the same application programming interface. ADflow is part of a wider aerodynamic shape optimization tool suite that is also available under an open-source license.

Robust Simulation of Coupled Reactions and Transport in Redox Flow Batteries Using Tailored Numerical Schemes

Redox flow batteries (RFBs) are attractive energy storage solutions for the grid, the simulation of which enables system and material optimization. In this article, we introduce tailored numerical schemes to model coupling between the transport of dissolved species, electrons, and fluid with redox reaction kinetics within RFBs in a robust way. The macro-scale transport of species (including advection, migration, and hydrodynamic dispersion) is coupled with the volume-averaged pore-scale processes of reaction kinetics, and mass transport, while the Poisson equation is used to model Donnan exclusion across ion exchange membranes that limit capacity fade due to the crossover of redox-active species. The governing equations are discretized using the finite volume method with a Newton-Raphson iteration scheme to resolve non-linearity. We introduce several numerical schemes to increase solver robustness, including reaction rate damping and logarithmic transform of concentration fields. Based on fixed-point iteration convergence criteria we also show that the mechanistic Marcus-Hush-Chidsey (MHC) redox kinetics model can tolerate larger time steps than the empirical Butler-Volmer (BV) kinetics model. The numerical schemes presented in this article can also find application in other electrochemical systems, including desalination devices, fuel cells, and electrodialysis.

Robust aerodynamic shape optimization—From a circle to an airfoil

Aerodynamic design optimization currently lacks robustness with respect to the starting design and requires trial and error in the flow solver and optimization algorithm settings to get a converged optimal design. We address need issue by developing ways to overcome robustness issues arising from shape parametrization, mesh deformation, and flow solver convergence. Our approach is demonstrated on the Aerodynamic Design Optimization Discussion Group (ADODG) airfoil optimization benchmarks to show the factors that dominate the robustness and efficiency. In the ADODG NACA 0012 benchmark, we address the additional issue of non-unique solutions. In the ADODG RAE 2822 case, we address solver failure due to shock waves, separation, and gradient accuracies due to the frozen turbulence model. Finally, we create a new, challenging aerodynamic shape optimization case that starts with a circle to test the robustness of our aerodynamic shape optimization framework. We use both fixed and adaptive parametrization methods to tackle this problem and show how we can exploit the advantages of adaptive parametrization methods to improve both robustness and efficiency. The combination of flow solver robustness, precision of gradient information, robust mesh deformation, and adaptive parametrization brings us closer to a “push-button” solution for airfoil design.

An accurate discontinuous Galerkin method for solving point-source Eikonal equation in 2-D heterogeneous anisotropic media

Accurate numerical computation of wave traveltimes in heterogeneous media is of major interest for a large range of applications in seismics, such as phase identification, data windowing, traveltime tomography and seismic imaging. A high level of precision is needed for traveltimes and their derivatives in applications which require quantities such as amplitude or take-off angle. Even more challenging is the anisotropic case, where the general Eikonal equation is a quartic in the derivatives of traveltimes. Despite their efficiency on Cartesian meshes, finite-difference solvers are inappropriate when dealing with unstructured meshes and irregular topographies. Moreover, reaching high orders of accuracy generally requires wide stencils and high additional computational load. To go beyond these limitations, we propose a discontinuous-finite-element-based strategy which has the following advantages: (1) the Hamiltonian formalism is general enough for handling the full anisotropic Eikonal equations; (2) the scheme is suitable for any desired high-order formulation or mixing of orders (p-adaptivity); (3) the solver is explicit whatever Hamiltonian is used (no need to find the roots of the quartic); (4) the use of unstructured meshes provides the flexibility for handling complex boundary geometries such as topographies (h-adaptivity) and radiation boundary conditions for mimicking an infinite medium. The point-source factorization principles are extended to this discontinuous Galerkin formulation. Extensive tests in smooth analytical media demonstrate the high accuracy of the method. Simulations in strongly heterogeneous media illustrate the solver robustness to realistic Earth-sciences-oriented applications.

Automatic generation of locally refined Cartesian meshes: Data management and algorithms

SummaryThe paper describes the main features of an automatic and three‐dimensional Cartesian mesher specifically designed for compressible inviscid/viscous flow solvers based on an immersed boundary technique. The development of a meshing tool able of conducting non‐isotropic cell refinements is a very tiresome task. The major difficulty is to imagine, at the pre‐design phase, a light but flexible data management, which minimizes the memory and CPUs' resources. In particular, the embedded geometry has to be detected by means of a fast and robust tagging procedure. Cells in proximity of the wall have to be refined in a proper way to adequately solve large flow gradients. Smooth variation of mesh density between differently refined zones has to be guaranteed to increase the flow solver robustness. A procedure to obtain accurate data on the geometry surfaces should be foresee. Here, a robust algorithm is developed to reconstruct a surface triangulation starting from the intersection points among volume cells and the geometry surfaces. The paper attempts to address all the above issues to help the readers in designing their own tools and suggesting them a way forward.

Application of high-order numerical schemes and Newton-Krylov method to two-phase drift-flux model

This paper concerns the application and solver robustness of the Newton-Krylov method in solving two-phase flow drift-flux model problems using high-order numerical schemes. In our previous studies, the Newton-Krylov method has been proven as a promising solver for two-phase flow drift-flux model problems. However, these studies were limited to use first-order numerical schemes only. Moreover, the previous approach to treating the drift-flux closure correlations was later revealed to cause deteriorated solver convergence performance, when the mesh was highly refined, and also when higher-order numerical schemes were employed. In this study, a second-order spatial discretization scheme that has been tested with two-fluid two-phase flow model was extended to solve drift-flux model problems. In order to improve solver robustness, and therefore efficiency, a new approach was proposed to treating the mean drift velocity of the gas phase as a primary nonlinear variable to the equation system. With this new approach, significant improvement in solver robustness was achieved. With highly refined mesh, the proposed treatment along with the Newton-Krylov solver were extensively tested with two-phase flow problems that cover a wide range of thermal-hydraulics conditions. Satisfactory convergence performances were observed for all test cases. Numerical verification was then performed in the form of mesh convergence studies, from which expected orders of accuracy were obtained for both the first-order and the second-order spatial discretization schemes. At last, the drift-flux model, along with numerical methods presented, were validated with three sets of flow boiling experiments that cover different flow channel geometries (round tube, rectangular tube, and rod bundle), and a wide range of test conditions (pressure, mass flux, wall heat flux, inlet subcooling and outlet void fraction).

A multidomain discretization of the Richards equation in layered soil

We consider the Richards equation on a domain that is decomposed into nonoverlapping layers, i.e., the decomposition has no cross points. We assume that the saturation and permeability functions are space-independent on each subdomain. Kirchhoff transformation of each subdomain problem separately then leads to a set of semilinear equations, which can each be solved efficiently using monotone multigrid. The transformed subdomain problems are coupled by nonlinear continuity and flux conditions. This nonlinear coupled problem can be solved using substructuring methods like the Dirichlet–Neumann or Robin iteration. We give several numerical examples showing the discretization error, the solver robustness under variations of the soil parameters, and a hydrological example with four soil layers and surface water.

Testing linear solvers for global gradient algorithm

Steady-state Water Distribution Network models compute pipe flows and nodal heads for assumed nodal demands, pipe hydraulic resistances, etc. The nonlinear mathematical problem is based on energy and mass conservation laws which is solved by using global linearization techniques, such as global gradient algorithm (GGA). The matrix of coefficients of the linear system inside GGA belongs to the class of sparse, symmetric and positive definite. Therefore a fast solver for the linear system is important in order to achieve the computational efficiency, especially when multiple runs are required. This work aims at testing three main strategies for the solution of linear systems inside GGA. The tests are performed on eight real networks by sampling nodal demands, considering the pressure-driven and demand-driven modelling to evaluate the robustness of solvers. The results show that there exists a robust specialized direct method which is superior to all the other alternatives. Furthermore, it is found that the number of times the linear system is solved inside the GGA does not depend on the specific solver, if a small regularization to the linear problem is applied, and that pressure-driven modelling requires a greater number which depends on the size and topology of the network and not only on the level of pressure deficiency.

A quantitative metric for robustness of nonlinear algebraic equation solvers

Abstract Practitioners in the area of dynamic simulation of technical systems report difficulties at times with steady-state initialization of models developed using general declarative modeling languages. These difficulties are analyzed in detail in this work and a rigorous approach to quantify robustness in the context of nonlinear algebraic equation systems is presented. This tool is then utilized in a study of six state of the art gradient-based iterative solvers on a set of industrial test problems. Finally, conclusions are drawn on the observed solver robustness in general, and it is argued whether the reported difficulties with steady-state initialization can be supported using the proposed quantitative metric.

Testing Computational Algorithms for Unsaturated Flow

The purpose of this work is to test different computational algorithms for unsaturated flow for accuracy and robustness by comparing computed results in a finite element program with analytical solutions. Because real-world problems are complex, testing codes for accuracy is often difficult. This is particularly true for flow in the vadose zone where Richards' equation is highly nonlinear. Recently, however, Tracy (Tracy WRRJ 2006) [1] (Tracy JHYD 2007) [2] has derived analytical solutions for a box-shaped flow region that is initially dry until water is applied to the top of the region. Two-dimensional and three-dimensional versions of these solutions for both steady-state and transient flow are available to be used in the testing process. Numerical precision and nonlinear solver robustness were investigated for varying degrees of nonlinearity by varying the Gardner parameter. As was increased, three ways of modeling relative hydraulic conductivity inside individual finite elements and two versions of the nonlinear solver were tested using three different ways to measure the error. The results of these tests are given in this paper.