Robust Solver Research Articles

Numerical simulation of carbon dioxide flooding in fractured reservoirs using generic projection-based embedded discrete fracture model

This paper, for the first time, integrates the generic projection-based embedded discrete fracture model (pEDFM) with the commercial reservoir simulator ECLIPSE for carbon dioxide (CO2) flooding in fractured reservoirs. The integrated method first obtains inter-grid connections and corresponding transmissibilities within the reservoir model based on the generic pEDFM. It then constructs an equivalent CO2 flooding ECLIPSE model to the original pEDFM reservoir model, thereby calculating the global equations of the compositional flow model to obtain distributions of pressure, saturation, component concentrations, and well performance data. We implemented three numerical examples to verify that the proposed method can achieve significantly higher computational accuracy compared to the widely used embedded discrete fracture model in both high and low permeability fracture scenarios, while also avoiding the difficulties associated with generating matching grids for complex fracture networks. Furthermore, the proposed integrated method uses the solver within ECLIPSE to solve the global equations, thus avoiding the high cost of developing a robust nonlinear solver for complex compositional model of CO2 flooding. This demonstrates the practicality of the method and its significant potential for subsequent application to various reservoir models.

3-D parallel anisotropic inversion of controlled-source electromagnetic data using nested tetrahedral grids

SUMMARY Geophysicists today face the challenge of quickly and reliably interpreting extensive controlled-source electromagnetic (CSEM) data sets to map subsurface conductivity structures within realistic geological environments. An ideal 3-D CSEM inversion algorithm using tetrahedral grids should be capable of distinguishing different resolution requirements between forward modelling and inversion grids, have an optimal parallel strategy that fully exploits the inherent independence of CSEM data sets while also possessing the capability to handle large-scale geo-electrical models, and incorporate conductivity anisotropy which should be a common characteristic in realistic subsurface environments. However, existing tools in the geo-electromagnetic community often fall short of these three demands. Addressing this gap, our study introduces a scalable and parallel anisotropic inversion technique for CSEM data, capitalizing on the potential of unstructured tetrahedral grids. We first apply the tetrahedral longest-edge bisection method to create a refined dense, heterogeneous forward modelling grid from a coarse inversion grid. This refinement, focused on areas around transmitters and receivers, is seamlessly integrated within the coarser inversion grid’s topology, enabling precise conductivity mapping and preserving electromagnetic response accuracy during model updates. We further innovate with a source-mesh double-level parallel strategy, utilizing the message passing interface technique for parallel handling of independent CSEM data sets and large-scale geo-electrical models. Externally, we dedicate a processor for inversion model updates employing the Limited-memory Broyden–Fletcher–Goldfarb–Shanno optimization algorithm and divide other processors into groups, each associated with specific transmitting sources and frequencies. Internally, in each group, we employ a domain-decomposition-based scalable and robust iterative solvers using the Auxiliary-Space Maxwell pre-conditioner to parallel quickly calculate the electromagnetic responses from its assigned source-frequency set. Additionally, recognizing the potential for electrical conductivity anisotropy in field data, we incorporate the case of vertical transverse isotropy. We validate the effectiveness of our method through examples, including an isotropic land model with undulating topography, an anisotropic marine model and a real-field data case. Results from both synthetic and field data inversions underscore our method’s significant advancements in efficiency and practicality, particularly in addressing large-scale 3-D CSEM data sets inversion challenges in realistic geological environments.

Revisiting Artifacts of Kohn-Sham Density Functionals for Biosimulation.

We revisit the problem of unphysical charge density delocalization/fractionalization induced by the self-interaction error of common approximate Kohn-Sham (KS) density functional theory functionals on simulation of small to medium-sized proteins in a vacuum. Aside from producing unphysical electron densities and total energies, the vanishing of the HOMO-LUMO gap associated with the unphysical charge delocalization leads to an unphysical low-energy spectrum and catastrophic failure of most popular solvers for the KS self-consistent field (SCF) problem. We apply a robust quasi-Newton SCF solver [ Phys. Chem. Chem. Phys. 2024, 26, 6557] to obtain solutions for some of these difficult cases. The anatomy of the charge delocalization is revealed by the natural deformation orbitals obtained from the density matrix difference between the Hartree-Fock and KS solutions; the charge delocalization not only can occur between charged fragments (such as in zwitterionic polypeptides) but also involves neutral fragments. The vanishing-gap phenomenon and troublesome SCF convergence are both attributed to the unphysical KS Fock operator eigenspectra of molecular fragments (e.g., amino acids or their side chains). Analysis of amino acid pairs suggests that the unphysical charge delocalization can be partially ameliorated by the use of some range-separated hybrid functionals but not by semilocal or standard hybrid functionals. Last, we demonstrate that solutions without the unphysical charge delocalization can be located even for semilocal KS functionals highly prone to such defects, but such solutions have non-Aufbau character and are unstable with respect to mixing of the non-overlapping "frontier" orbitals. Caution should be exercised when unexpectedly small (or vanishing) HOMO-LUMO gaps and atypical SCF convergence patterns (e.g., oscillatory) are observed in KS DFT simulations in any context (bio or otherwise).

A robust and well-balanced finite volume solver for investigating the effects of tides on water renewal timescale in the Nador lagoon, Morocco

Coastal lagoons, particularly those not well-connected to the sea, are highly susceptible to continuous water retention, adversely affecting water quality and ecosystems. Understanding and quantifying the processes governing water renewal in such semi-enclosed regions is crucial. Using the principles of constituent-oriented age and residence time theory (www.climate.be/cart), we calculated timescales to explain the water renewal process in semi-enclosed areas. To improve our understanding and define the dynamics of the Nador lagoon, we use two distinct time scales: residence time and exposure time. Residence time indicates the length of time a parcel of water leaves the region of interest for the first time, while exposure time represents the overall length of time a parcel of water remains in that region, including any subsequent inflows. The modeling system adopts an Eulerian approach, including two interlinked model components: the shallow-water equations incorporating bottom friction, eddy viscosity, wind shear stresses, momentum diffusion, and Coriolis forces for modeling hydrodynamics, and a transport-diffusion equation utilized to simulate the advection and diffusion of the passive tracer concentration. The entire system is integrated into a high-order finite volume solver that employs unstructured meshes, upwind numerical fluxes, and slope limiters to achieve precise resolution of steep bathymetric gradients that may arise in the approximate solution. The scheme is non-oscillatory and preserves the conservation properties. In this research, we explore various situations related to the connection between the Mediterranean Sea and the Nador lagoon. Specifically, we suggest novel entry points that have the potential to significantly decrease the water residence within the lagoon.

URANOS-2.0: Improved performance, enhanced portability, and model extension towards exascale computing of high-speed engineering flows

We present URANOS-2.0, the second major release of our massively parallel, GPU-accelerated solver for compressible wall flow applications. This latest version represents a significant leap forward in our initial tool, which was launched in 2023 (De Vanna et al. [1]), and has been specifically optimized to take full advantage of the opportunities offered by the cutting-edge pre-exascale architectures available within the EuroHPC JU. In particular, URANOS-2.0 emphasizes portability and compatibility improvements with the two top-ranked supercomputing architectures in Europe: LUMI and Leonardo. These systems utilize different GPU architectures, AMD and NVIDIA, respectively, which necessitates extensive efforts to ensure seamless usability across their distinct structures. In pursuit of this objective, the current release adheres to the OpenACC standard. This choice not only facilitates efficient utilization of the full potential inherent in these extensive GPU-based architectures but also upholds the principles of vendor neutrality, a distinctive characteristic of URANOS solvers in the CFD solvers' panorama. However, the URANOS-2.0 version goes beyond the goals of improving usability and portability; it introduces performance enhancements and restructures the most demanding computational kernels. This translates into a 2× speedup over the same architecture. In addition to its enhanced single-GPU performance, the present solver release demonstrates very good scalability in multi-GPU environments. URANOS-2.0, in fact, achieves strong scaling efficiencies of over 80% across 64 compute nodes (256 GPUs) for both LUMI and Leonardo. Furthermore, its weak scaling efficiencies reach approximately 95% and 90% on LUMI and Leonardo, respectively, when up to 256 nodes (1024 GPUs) are considered. These significant performance advancements position URANOS-2.0 as a state-of-the-art supercomputing platform tailored for compressible wall turbulence applications, establishing the solver as an integrated tool for various aerospace and energy engineering applications, which can span from direct numerical simulations, wall-resolved large eddy simulations, up to most recent wall-modeled large eddy simulations. Program summaryProgram title: Unsteady Robust All-around Navier-StOkes Solver (URANOS)CPC Library link to program files:https://doi.org/10.17632/pw5hshn9k6.2Developer's repository link:https://github.com/uranos-gpu/uranos-gpu, https://github.com/uranos-gpu/uranos-gpu/tree/v2.0Licensing provisions: BSD License 2.0Programming language: Modern Fortran, OpenACC, MPINature of problem: Solving the compressible Navier-Stokes equations in a three-dimensional Cartesian framework.Solution method: Convective terms are treated with high-resolution shock-capturing schemes. The system dynamics is advanced in time with a three-stage Runge-Kutta method. Parallelization adopts MPI+OpenACC.

Turbulence driven goal-oriented anisotropic mesh adaptation for RANS simulations in aerodynamics

Anisotropic mesh adaptation is starting to play a role of great importance within the CFD community. In fact, such a tool allows to accurately capture the physics of the problem in an automatic way, by reducing at the same time the computational resources, in contrast to the best-practice meshing guidelines, which require an a priori knowledge of the flow solution and, for this reason, are approximate and error-prone. In the present work, we consider anisotropic mesh adaptation for steady flows to get mesh-independent solutions with a small number of vertices, that is, we achieve the so-called early capturing property. Two classes of methods exist: feature-based mesh adaptation where the minimization of the interpolation error associated with sensor field is targeted and goal-oriented mesh adaptation where the minimization of the discretization error associated with a goal functional is carried out, such as the lift or the drag coefficients. The goal-oriented mesh adaptation requires the adjoint state computation, thus a crucial part of the treatment deals with this topic. The novelty of the present approach is the taking into account of the turbulent primal and adjoint variables, and the corresponding fluxes (including turbulence sources) inside the mesh adaptation process and the adjoint computation, as up to now only the inviscid and the viscous contributions have been considered in this kind of processes. We show how such additions are extremely beneficial for goal-oriented mesh adaptation, as the resulting meshes are able to capture the physics of the problem with a higher accuracy, in contrast to the exclusively inviscid- or viscous-driven adaptation. Eventually, as the mesh adaptation process tends to produce highly unstructured anisotropic meshes, a relevant part of the present work is devoted to the description of a robust solver capable of solving the RANS equations on such meshes.

A mixed nonlocal finite element model for thermo‐poro‐elasto‐plastic simulation of porous media with multiphase fluid flow

SummaryA new mixed nonlocal finite element framework is developed for nonlinear thermo‐poro‐elasto‐plastic simulation of porous media with multiphase pore fluid flow and thermal coupling. The solid‐fluid interaction is accounted for using the mixture theory of Biot based on the volume fractions concept. Different sources of nolinearities arising from the multiphase fluid flow effects, advective‐diffusive heat transfer, inelastic deformation, fluid flux injection induced mechanical tractions, solid skeleton deformation permeability dependence, and temperature dependent viscosity are included in developing a robust numerical solver for the targeted coupled multiphysics problem. To address the effect of microstructure in inelastic localized deformation behaviour with dilational softening, a nonlocal plasticity model is proposed based on a characteristic length scale which rectifies the non‐physical pathological mesh dependence problem encountered in conventional plasticity. The accuracy and strength of the developed model is shown with comparing the obtained numerical results of a benchmark bilateral compression test with existing published data in the literature. To show the versatility and robustness of the developed computational framework in modelling the geomechanics of real‐case engineering practices, large scale thermo‐hydro‐mechanical (THM) subsurface stimulation processes with applications in enhanced oil recovery (EOR) are effectively simulated and the targeted enhanced recovery and performances are demonstrated. The current formulation does not include phase transformation modelling capability, and therefore, the developed models may not be applicable for simulation of the engineering processes that involve phase change behaviour (e.g., steam injection).

Recycling Krylov subspaces for efficient partitioned solution of aerostructural adjoint systems

Robust and efficient solvers for coupled-adjoint linear systems are crucial to successful aerostructural optimization. Monolithic and partitioned strategies can be applied. The monolithic approach is expected to offer better robustness and efficiency for strong fluid-structure interactions. However, it requires a higher implementation cost and convergence may depend on appropriate scaling and initialization strategies. On the other hand, the modularity of the partitioned method enables a straightforward implementation while its convergence may require relaxation. In addition, a partitioned solver often leads to a higher number of iterations to get the same level of convergence as the monolithic one.The objective of this paper is to accelerate the fluid-structure coupled-adjoint partitioned solver by considering techniques borrowed from approximate invariant subspace recycling strategies adapted to sequences of linear systems with varying right-hand sides. Indeed, in a partitioned framework, the structural source term attached to the fluid block of equations affects the right-hand side with the nice property of quickly converging to a constant value. For the fluid block of equations, we also consider deflation of approximate eigenvectors in conjunction with advanced inner-outer Krylov solvers.To demonstrate the effectiveness of our proposed recycling strategy, we compute the coupled derivatives for two emblematic problems in transonic viscous flow. The first one is an aeroelastic configuration of the ONERA-M6 fixed wing. For this exercise, the fluid grid was coupled to a structural model specifically designed to exhibit a high flexibility. The second one is a wing-body aeroelastic configuration of the Common Research Model. All computations are performed using a fully linearized one-equation Spalart-Allmaras turbulence model. Numerical simulations show up to 39% reduction in matrix-vector products for GCRO-DR and up to 25% for the nested FGCRO-DR solver.

DESIGN AND ANALYSIS OF THE PISTON USING COMPOSITE MATERIALS LIKE AA7275 STRUCTURAL AND THERMAL ANALYSIS DONE IN ANSYS SOFTWARE

This project focuses on the design and analysis of a piston utilizing composite materials, particularly AA7275, to enhance structural integrity and thermal efficiency. The study employs ANSYS software to conduct comprehensive structural and thermal analyses. The structural analysis evaluates the mechanical behavior, stress distribution, and deformation characteristics of the piston under varying operating conditions. Meanwhile, the thermal analysis investigates the heat transfer mechanisms, temperature distribution, and thermal stress accumulation within the piston. By integrating composite materials and utilizing advanced simulation techniques, this project aims to optimize the performance and reliability of the piston, contributing to advancements in engine technology. ANSYS software facilitates Multiphysics simulations, allowing engineers to study the interactions between different physical phenomena such as structural mechanics, fluid dynamics, electromagnetics, and thermal effects. Its robust solvers and algorithms enable accurate prediction of system behavior under various operating conditions, helping engineers optimize performance, reliability, and efficiency. Moreover, ANSYS provides a range of specialized modules tailored to specific industries and applications, including automotive, aerospace, electronics, healthcare, and renewable energy. These modules incorporate domain-specific features and workflows, empowering engineers to address industry-specific challenges and design requirements effectively. Recent advancements in ANSYS software include enhanced computational capabilities, improved user interfaces, and integration with emerging technologies such as artificial intelligence and additive manufacturing. These developments enable engineers to accelerate the design iteration process, reduce time-to-market, and achieve superior product performance. Furthermore, ANSYS supports collaborative engineering workflows through its integration with product lifecycle management (PLM) and computer-aided design (CAD) software. This integration facilitates seamless data exchange between design, simulation, and manufacturing teams, ensuring consistency and accuracy throughout the product development lifecycle. KEY WORDS : AA7275; Piston; Structural integrity and Thermal Efficiency; Structural analysis; ANSYS; CAD.

Compositional Simulation for Carbon Storage in Porous Media Using an Electrolyte Association Equation of State

Summary We present a new algorithm based on automatic differentiation that enables precise computation of the derivatives of the Z-factor, facilitating the utilization of Newton’s method or coupling with a robust flow solver. Leveraging a free open-source code [MATLAB Reservoir Simulation Toolbox (MRST)], we develop an electrolyte cubic plus association (e-CPA) equation of state (EoS) model to accurately represent the injection of carbon dioxide (CO2) in brine. By integrating flow and thermodynamics, we construct an advanced compositional simulator using MRST’s object-oriented, automatic differentiation framework and the newly developed e-CPA EoS model. This simulator offers flexibility through both overall-composition and natural-variable formulations, achieved by selecting different primary variables. The Péneloux volume translation technique is employed to modify the EoS model’s volume, ensuring accurate density calculation for the mixture. Additionally, we introduce a viscosity model, e-CPA-FV, which accurately predicts the viscosity of carbon capture and storage (CCS) fluids, surpassing the accuracy of the traditional Lohrenz-Bray-Clark (LBC) model. Our simulator demonstrates superior performance in predicting CO2-brine systems compared with the standard formulation based on the Peng-Robinson (PR) EoS and can handle brine with various salts. The self-contained source code necessary to reproduce all examples is available on the open-access Zenodo digital repository (doi: 10.5281/zenodo.10691505).

Robust and accurate Roe-type Riemann solver with compact stencil: Rotated-RoeM scheme

In this study, we aim to develop a robust and accurate Riemann solver to resolve strong shock waves with compact stencil information by applying the rotated hybrid concept to the RoeM scheme. The original RoeM scheme offers good robustness from the perspective of avoiding shock instability, also called the carbuncle phenomenon, in hypersonic flows by introducing Mach-number-based functions. These functions rely on multi-dimensional shock discontinuity sensing terms to capture grid-aligned shock waves, thus necessitating sharing of values of cells adjacent to each vertex. Consequently, constructing a data structure that can handle such a wide stencil is necessary, particularly for implementations on unstructured grids. Moreover, communication of this stencil information on parallel computation can increase computational costs.To overcome these complexities, the rotated hybrid concept is applied to the RoeM scheme, to replace the multi-dimensional terms while still resolving the grid-aligned shock instability. In addition, a Mach-number-based weighting function is introduced to distinguish between the flow around the boundary layer and the onset of shock instability. The damping characteristic of the proposed scheme is compared to that of the original RoeM scheme by conducting a linear perturbation analysis. The findings reveal that the proposed scheme effectively controls spurious oscillations across strong shock waves. Numerous test cases are presented to demonstrate the robustness and accuracy of the proposed approach. Notably, the proposed scheme successfully resolves thermal boundary layers and accurately predicts the heat transfer rate around stagnation points in hypersonic flows.

Color-gradient-based phase-field equationfor multiphase flow.

In this paper, the underlying problem with the color-gradient (CG) method in handling density-contrast fluids is explored. It is shown that the CG method is not fluid invariant. Based on nondimensionalizing the CG method, a phase-field interface-capturing model is proposed which tackles the difficulty of handling density-contrast fluids. The proposed formulation is developed for incompressible, immiscible two-fluid flows without phase-change phenomena, and a solver based on the lattice Boltzmann method is proposed. Coupled with an available robust hydrodynamic solver, a binary fluid flow package that handles fluid flows with high density and viscosity contrasts is presented. The macroscopic and lattice Boltzmann equivalents of the formulation, which make the physical interpretation of it easier, are presented. In contrast to existing color-gradient models where the interface-capturing equationsare coupled with the hydrodynamic ones and include the surface tension forces, the proposed formulation is in the same spirit as the other phase-field models such as the Cahn-Hilliard and the Allen-Cahn equationsand is solely employed to capture the interface advected due to a flow velocity. As such, similarly to other phase-field models, a so-called mobility parameter comes into play. In contrast, the mobility is not related to the density field but a constant coefficient. This leads to a formulation that avoids individual speed of sound for the different fluids. On the lattice Boltzmann solver side, two separate distribution functions are adopted to solve the formulation, and another one is employed to solve the Navier-Stokes equations, yielding a total of three equations. Two series of numerical tests are conducted to validate the accuracy and stability of the model, where we compare simulated results with available analytical and numerical solutions, and good agreement is observed. In the first set the interfacial evolution equationsare assessed, while in the second set the hydrodynamic effects are taken into account.

A Positivity-Preserving and Robust Fast Solver for Time-Fractional Convection–Diffusion Problems

A Positivity-Preserving and Robust Fast Solver for Time-Fractional Convection–Diffusion Problems

A Robust and Scalable Multigrid Solver for 3-D Low-Frequency Electromagnetic Diffusion Problems

A Robust and Scalable Multigrid Solver for 3-D Low-Frequency Electromagnetic Diffusion Problems

Hp-Robust multigrid solver on locally refined meshes for FEM discretizations of symmetric elliptic PDEs

In this work, we formulate and analyze a geometric multigrid method for the iterative solution of the discrete systems arising from the finite element discretization of symmetric second-order linear elliptic diffusion problems. We show that the iterative solver contracts the algebraic error robustly with respect to the polynomial degree p ≥ 1 and the (local) mesh size h. We further prove that the built-in algebraic error estimator which comes with the solver is hp-robustly equivalent to the algebraic error. The application of the solver within the framework of adaptive finite element methods with quasi-optimal computational cost is outlined. Numerical experiments confirm the theoretical findings.

Generalizing reduction‐based algebraic multigrid

AbstractAlgebraic multigrid (AMG) methods are often robust and effective solvers for solving the large and sparse linear systems that arise from discretized PDEs and other problems, relying on heuristic graph algorithms to achieve their performance. Reduction‐based AMG (AMGr) algorithms attempt to formalize these heuristics by providing two‐level convergence bounds that depend concretely on properties of the partitioning of the given matrix into its fine‐ and coarse‐grid degrees of freedom. MacLachlan and Saad (SISC 2007) proved that the AMGr method yields provably robust two‐level convergence for symmetric and positive‐definite matrices that are diagonally dominant, with a convergence factor bounded as a function of a coarsening parameter. However, when applying AMGr algorithms to matrices that are not diagonally dominant, not only do the convergence factor bounds not hold, but measured performance is notably degraded. Here, we present modifications to the classical AMGr algorithm that improve its performance on matrices that are not diagonally dominant, making use of strength of connection, sparse approximate inverse (SPAI) techniques, and interpolation truncation and rescaling, to improve robustness while maintaining control of the algorithmic costs. We present numerical results demonstrating the robustness of this approach for both classical isotropic diffusion problems and for non‐diagonally dominant systems coming from anisotropic diffusion.

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.

Hotspot resolution in cloud computing: A Γ-robust knapsack approach for virtual machine migration

Cloud providers offer virtual machines (VMs) located on physical machines (PMs) to meet the increasing demand for computational services. When the instantaneous utilized capacities of VMs exceed a PM's threshold, hotspots form and degrade VM performance. To resolve hotspots, some VMs must be live migrated to other PMs, but selecting which VMs is challenging as their utilized capacities change continuously. We propose a Predicted Mixed Integer Linear Programming (MILP) Robust Solver (PMRS) that applies Γ-robustness theory to ensure PMs are hotspot-safe with a desired probability. PMRS uses a “prediction + optimization” framework that first predicts VMs' future behaviors and then formulates the problem as a Γ-robust knapsack problem (Γ-RKP) solvable with a novel MILP model. Experiments with real-trace and synthetic data demonstrate PMRS's effectiveness. Moreover, we apply PMRS in a real production environment in Huawei Cloud, and observe significant benefits in resolving existing hotspots and 94%+ potential future hotspots with minimal migration cost.

Comparison of Lattice Boltzmann and Boundary Element Methods for Stokes and Visco-Inertial Flow in a Two-Dimensional Porous Medium

In porous media, limitations imposed by macroscale laws can be avoided with a dual-scale model, in which the pore-scale phenomena of interest are modelled directly over a large number of realisations. Such a model requires a robust, accurate and efficient pore-scale solver. We compare the boundary element method (BEM) and two variants of the lattice Boltzmann method (LBM) as pore-scale solvers of 2D incompressible flow. The methods are run on a number of test cases and the performance of each simulation is assessed according to the mean velocity error and the computational runtime. Both the porous geometry (porosity, shape and complexity), and the Reynolds number (from Stokes to visco-inertial flow) are varied between the test cases. We find that, for Stokes flow, BEM provides the most efficient and accurate solution in simple geometries (with small boundary length) or when a large runtime is practical. In all other situations we consider, one of the variants of LBM performs best. We furthermore demonstrate that these findings also apply in a dual-scale model of Stokes flow through a locally-periodic medium.

Implementation of Asynchronous Distributed Gauss-Newton Optimization Algorithms for Uncertainty Quantification by Conditioning to Production Data

Summary Previous implementation of the distributed Gauss-Newton (DGN) optimization algorithm ran multiple optimization threads in parallel, employing a synchronous running mode (S-DGN). As a result, it waits for all simulations submitted in each iteration to complete, which may significantly degrade performance because a few simulations may run much longer than others, especially for time-consuming real-field cases. To overcome this limitation and thus improve the DGN optimizer’s execution, we propose two asynchronous DGN (A-DGN) optimization algorithms in this paper. The two A-DGN optimization algorithms are (1) the local-search algorithm (A-DGN-LS) to locate multiple maximum a-posteriori (MAP) estimates and (2) the integrated global-search algorithm with the randomized maximum likelihood (RML) method (A-DGN + RML) to generate hundreds of RML samples in parallel for uncertainty quantification. We propose using batch together with a checking time interval to control the optimization process. The A-DGN optimizers check the status of all running simulations after every checking time frame. The iteration index of each optimization thread is updated dynamically according to its simulation status. Thus, different optimization threads may have different iteration indices in the same batch. A new simulation case is proposed immediately once the simulation of an optimization thread is completed, without waiting for the completion of other simulations. We modified the training data set updating algorithm using each thread’s dynamically updated iteration index to implement the asynchronous running mode. We apply the modified QR decomposition method to estimate the sensitivity matrix at the best solution of each optimization thread by linear interpolation of all or a subset of the training data to avoid the issue of solving a linear system with a singular matrix because of insufficient training data points in early batches. A new simulation case (or search point) is generated by solving the Gauss-Newton (GN) trust-region subproblem (GNTRS) using the estimated sensitivity matrix. We developed a more efficient and robust GNTRS solver using eigenvalue decomposition (EVD). The proposed A-DGN optimization methods are tested and validated on a 2D analytical toy problem and a synthetic history-matching problem and then applied to a real-field deepwater reservoir model. Numerical tests confirm that the proposed A-DGN optimization methods can converge to solutions with matching quality comparable to those obtained by the S-DGN optimizers, saving on the time required for the optimizer to converge by a factor ranging from 1.3 to 2 when compared to the S-DGN optimizer depending on the problem. The new A-DGN optimization algorithms improve efficiency and robustness in solving history-matching or inversion problems, especially for uncertainty quantification of subsurface model parameters and production forecasts of real-field reservoirs by conditioning production data.