Anderson acceleration (AA) has a long history of use and a strong recent interest due to its potential ability to dramatically improve the linear convergence of the fixed-point iteration. Most authors are simply using and analyzing the stationary version of Anderson acceleration (sAA) with a constant damping factor or without damping. Little attention has been paid to nonstationary algorithms. However, damping can be useful and is sometimes crucial for simulations in which the underlying fixed-point operator is not globally contractive. The role of this damping factor has not been fully understood. In the present work, we consider the non-stationary Anderson acceleration algorithm with optimized damping (AAoptD) in each iteration to further speed up linear and nonlinear iterations by applying one extra inexpensive optimization. We analyze the convergence rate this procedure and develop an efficient and inexpensive implementation scheme. We show by extensive numerical experiments that the proposed non-stationary Anderson acceleration with optimized damping procedure often converges much faster than stationary AA with constant damping, adaptive damping or without damping, especially in the cases larger window sizes are needed. We also observe that simple strategies like using constant damping factors and adaptive damping factors, sometimes, work very well for some problems while sometimes they are even worse than AA without damping. Our proposed method is usually more robust than AA with constant damping and adaptive damping. Moreover, we also observed from our numerical results that damping can be good, but choosing the wrong damping factors may slow down the convergence rate. Theoretical analysis of the effects of damping factors are needed and important.

@en

Background: Individuals who were formerly incarcerated have high tuberculosis incidence, but are generally not considered among the risk groups eligible for tuberculosis prevention. We investigated the potential health impact and cost-effectiveness of Mycobacterium tuberculosis infection screening and tuberculosis preventive treatment (TPT) for individuals who were formerly incarcerated in Brazil. Methods: Using published evidence for Brazil, we constructed a Markov state transition model estimating tuberculosis-related health outcomes and costs among individuals who were formerly incarcerated, by simulating transitions between health states over time. The analysis compared tuberculosis infection screening and TPT, to no screening, considering a combination of M tuberculosis infection tests and TPT regimens. We quantified health effects as reductions in tuberculosis cases, tuberculosis deaths, and disability-adjusted life-years (DALYs). We assessed costs from a tuberculosis programme perspective. We report intervention cost-effectiveness as the incremental costs per DALY averted, and tested how results changed across subgroups of the target population. Findings: Compared with no intervention, an intervention incorporating tuberculin skin testing and treatment with 3 months of isoniazid and rifapentine would avert 31 (95% uncertainty interval 14–56) lifetime tuberculosis cases and 4·1 (1·4–5·8) lifetime tuberculosis deaths per 1000 individuals, and cost US$242 per DALY averted. All test and regimen combinations were cost-effective compared with no screening. Younger age, longer incarceration, and more recent prison release were each associated with significantly greater health benefits and more favourable cost-effectiveness ratios, although the intervention was cost-effective for all subgroups examined. Interpretation: M tuberculosis infection screening and TPT for individuals who were formerly incarcerated appears cost-effective, and would provide valuable health gains. Funding: National Institutes of Health. Translation: For the Portuguese translation of the abstract see Supplementary Materials section.

@en

In this paper, two new iterative methods for solving generalized absolute value equations (GAVEs) are proposed and investigated using the single-step iteration (SSI) approach. The proposed iterative methods are Picard-SSI and nonlinear SSI-like methods. In the implementation of the Picard-SSI method, we have used the SSI method as an inner solver. The convergence of the proposed method for solving GAVE is analyzed under reasonable constraints. Several numerical examples are given to illustrate the efficiency and implementation of the proposed methods.@en

We propose a matrix-free parallel two-level deflation method combined with the Complex Shifted Laplacian Preconditioner (CSLP) for two-dimensional heterogeneous Helmholtz problems encountered in seismic exploration, antennas, and medical imaging. These problems pose challenges in terms of accuracy and convergence due to scalability issues with numerical solvers. Motivated by the limitations imposed by excessive computational time and memory constraints when employing a sequential solver with constructed matrices, we parallelize the two-level deflation method without constructing any matrices. Our approach utilizes preconditioned Krylov subspace methods and approximates the CSLP preconditioner with a parallel geometric multigrid V-cycle. For the two-level deflation, standard inter-grid deflation vectors and further high-order deflation vectors are considered. As another main contribution, the matrix-free Galerkin coarsening approach and a novel re-discretization scheme as well as high-order finite-difference schemes on the coarse grid are studied to obtain wavenumber-independent convergence. The optimal settings for an efficient coarse-grid problem solver are investigated. Numerical experiments of model problems show that the wavenumber independence has been obtained for medium wavenumbers. The matrix-free parallel framework shows satisfactory weak and strong parallel scalability.

@en

We present a matrix-free parallel iterative solver for the Helmholtz equation related to applications in seismic problems and study its parallel performance. We apply Krylov subspace methods, GMRES, Bi-CGSTAB and IDR(s), to solve the linear system obtained from a second-order finite difference discretization. The Complex Shifted Laplace Preconditioner (CSLP) is employed to improve the convergence of Krylov solvers. The preconditioner is approximately inverted by multigrid iterations. For parallel computing, the global domain is partitioned blockwise. The standard MPI library is employed for data communication. The matrix-vector multiplication and preconditioning operator are implemented in a matrix-free way instead of constructing large, memory-consuming coefficient matrices. These adjustments lead to direct improvements in terms of memory consumption. Numerical experiments of model problems show that the matrix-free parallel solution method has satisfactory parallel performance and weak scalability. It allows us to solve larger problems in parallel to obtain more accurate numerical solutions.@en

The Helmholtz equation is related to seismic exploration, sonar, antennas, and medical imaging applications. It is one of the most challenging problems to solve in terms of accuracy and convergence due to the scalability issues of the numerical solvers. For 3D large-scale applications, high-performance parallel solvers are also needed. In this paper, a matrix-free parallel iterative solver is presented for the three-dimensional (3D) heterogeneous Helmholtz equation. We consider the preconditioned Krylov subspace methods for solving the linear system obtained from finite-difference discretization. The Complex Shifted Laplace Preconditioner (CSLP) is employed since it results in a linear increase in the number of iterations as a function of the wavenumber. The preconditioner is approximately inverted using one parallel 3D multigrid cycle. For parallel computing, the global domain is partitioned blockwise. The matrix-vector multiplication and preconditioning operator are implemented in a matrix-free way instead of constructing large, memory-consuming coefficient matrices. Numerical experiments of 3D model problems demonstrate the robustness and outstanding strong scaling of our matrix-free parallel solution method. Moreover, the weak parallel scalability indicates our approach is suitable for realistic 3D heterogeneous Helmholtz problems with minimized pollution error.@en

This study explores the suitability of quasi-static pore-network modeling for simulating the transport of hydrogen in networks with box-shaped pores and square cylinder throats. The dynamic pore-network modeling results are compared with quasi-static pore-network modeling, and a good agreement is observed when the simulations reach steady-state, for a capillary number of N_c≤10⁻⁷. This finding suggests that the quasi-static approach can be used as a reliable and efficient method for studying hydrogen transport in similar networks.

@en

Modeling of fluid flow in porous media is a pillar in geoscience applications. Previous studies have revealed that heterogeneity and fracture distribution have considerable influence on fluid flow. In this work, a numerical investigation of two-phase flow in heterogeneous fractured reservoir is presented. First, the discrete fracture model is implemented based on a hybrid-dimensional modeling approach, and an equivalent continuum approach is integrated in the model to reduce computational cost. A multilevel adaptive strategy is devised to improve the numerical robustness and efficiency. It allows up to 4-levels adaption, where the adaptive factors can be modified flexibly. Then, numerical tests are conducted to verify the the proposed method and to evaluate its performance. Different adaptive strategies with 3-levels, 4-levels and fixed time schemes are analyzed to evaluate the computational cost and convergence history. These evaluations demonstrate the merits of this method compared to the classical method. Later, the heterogeneity in permeability field, as well as initial saturation, is modeled in a layer model, where the effect of layer angle and permeability on fluid flow is investigated. A porous medium containing multiple length fractures with different distributions is simulated. The fine-scale fractures are upscaled based on the equivalent approach, while the large-scale fractures are retained. The conductivity of the rock matrix is enhanced by the upscaled fine-scale fractures. The difference of hydraulic property between homogeneous and heterogeneous situations is analyzed. It reveals that the heterogeneity may influence fluid flow and production, while these impacts are also related to fracture distribution and permeability.

@en

In this book we discuss several numerical methods for solving ordinary differential equations. We emphasize the aspects that play an important role in practical problems. We confine ourselves to ordinary differential equations with the exception of the last chapter in which we discuss the heat equation, a parabolic partial differential equation. The techniques discussed in the intro-ductory chapters, for instance interpolation, numerical quadrature and the solution to nonlinear equations, may also be used outside the context of differential equations. They have been in-cluded to make the book self-contained as far as the numerical aspects are concerned. Chapters, sections and exercises marked with a * are not part of the Delft Institutional Package. The numerical examples in this book were implemented in Matlab, but also Python or any other programming language could be used. A list of references to background knowledge and related literature can be found at the end of this book. Extra information about this course can be found at http://NMODE.ewi.tudelft.nl, among which old exams, answers to the exercises, and a link to an online education platform. We thank Matthias Moller for his thorough reading of the draft of this book and his helpful suggestions.@en

In this paper, we consider a block Jacobi preconditioner and various deflation techniques applied in the Deflated Preconditioned Conjugate Gradient (DPCG) method for solving a sparse system of linear equations derived from a statistical linear mixed model that analyses simultaneously phenotypic and pedigree information of genotyped and ungenotyped animals with Single Polymorphism Nucleotide genotypes of genotyped animals. In livestock production systems, evaluating the genetic merit of the animals through such a model is a key process to ensure an improvement of animals for some characteristics of interest at each generation. First, we propose to define the deflation vectors using a subdomain deflation approach that considers some biological properties of the genotypes. Using simulated data, this approach reduces the number of iterations by up to 87% in comparison to a Preconditioned Conjugate Gradient method with a Jacobi preconditioner. Furthermore, compared to a DPCG method with same number of subdomains but defined randomly, this approach reduces the number of iterations by up to 20% for the same computational costs of one DPCG iteration. The properties of the resulting systems show that this approach annihilates the largest eigenvalues of the preconditioned coefficient matrix. Second, we propose the use of solution vectors of 12 systems of equations that include between 0.25% and 3% less data, as deflation vectors. For reducing the computational costs, we also consider a Proper Orthogonal Decomposition-reduced set of these 12 vectors. The properties of the resulting systems show that this recycling information approach annihilates the smallest eigenvalues of the preconditioned coefficient matrix, and results in a reduction of up to 39% in comparison to the PCG method. Finally, based on our experiment, the combination of the subdomain deflation approach relying on biological properties and of the POD-based approach to recycle previous solution vectors, for defining the deflation vectors, results in annihilating both the smallest and largest eigenvalues, and in a reduction of up to 88 % of the number of iterations in comparison to the PCG method.

@en

The accuracy, stability and computational efficiency of numerical methods on central processing units (CPUs) for the depth-averaged shallow water equations were well covered in the literature. A large number of these methods were already developed and compared. However, on graphics processing units (GPUs), such comparisons are relatively scarce. In this paper, we present the results of comparing two time-integration methods for the shallow water equations on structured grids. An explicit and a semi-implicit time integration method were considered. For the semi-implicit method, the performance of several iterative solvers was compared. The implementation of the semi-implicit method on a GPU in this study was a novel approach for the shallow water equations. This also holds for the repeated red black (RRB) solver that was found to be very efficient on a GPU. Additionally, the results of both methods were compared with several CPU-based software systems for the shallow water flows on structured grids. On a GPU, the simulations were 25 to 75 times faster than on a CPU. Theory predicts an explicit method to be best suited for a GPU due to the higher level of inherent parallelism. It was found that both the explicit and the semi-implicit methods ran efficiently on a GPU. For very shallow applications, the explicit method was preferred because the stability condition on the time step was not very restrictive. However, for deep water applications, we expect the semi-implicit method to be preferred.@en

Although Anderson acceleration AA(m) has been widely used to speed up nonlinear solvers, most authors are simply using and studying the stationary version of Anderson acceleration. The behavior and full potential of the non-stationary version of Anderson acceleration methods remain an open question. Motivated by the hybrid linear solver GMRESR (GMRES Recursive), we recently proposed a set of non-stationary Anderson acceleration algorithms with dynamic window sizes AA(m,AA(n)) for solving both linear and nonlinear problems. Significant gains are observed for our proposed algorithms but these gains are not well understood. In the present work, we first consider the case of using AA(m,AA(1)) for accelerating linear fixed-point iteration and derive the polynomial residual update formulas for non-stationary AA(m,AA(1)). Like stationary AA(m), we find that AA(m,AA(1)) with general initial guesses is also a multi-Krylov method and possesses a memory effect. However, AA(m,AA(1)) has higher order degree of polynomials and a stronger memory effect than that of AA(m) at the k-th iteration, which might explain the better performance of AA(m,AA(1)) compared to AA(m) as observed in our numerical experiments. Moreover, we further study the influence of initial guess on the asymptotic convergence factor of AA(1, AA(1)). We show a scaling invariance property of the initial guess x for the AA(1,AA(1)) method in the linear case. Then, we study the root-linear asymptotic convergence factor under scaling of the initial guess and we explicitly indicate the dependence of root-linear asymptotic convergence factors on the initial guess. Lastly, we numerically examine the influence of the initial guess on the asymptotic convergence factor of AA(m) and AA(m,AA(n)) for both linear and nonlinear problems.

@en

An efficient compositional framework is developed for simulation of CO ₂ storage in saline aquifers during a full-cycle injection, migration and post-migration processes. Essential trapping mechanisms, including structural, dissolution, and residual trapping, which operate at different time scales, are accurately captured in the presented unified framework. In particular, a parameterization method is proposed to efficiently describe the relevant physical processes. The proposed framework is validated by comparing the dynamics of gravity-induced convective transport with that reported in the literature. Results show good agreement for both the characteristics of descending fingers and the associated dissolution rate. The developed simulator is then applied to study the FluidFlower benchmark model. An experimental setup with heterogeneous geological layers is discretized into a two-dimensional computational domain where numerical simulation is performed. Impacts of hysteresis and the diffusion of CO ₂ in liquid phase on the migration and trapping of CO ₂ plume are investigated. Inclusion of the hysteresis effect does not affect plume migration in this benchmark model, whereas diffusion plays an important role in promoting convective mixing. This work casts a promising approach to predict the migration of the CO ₂ plume, and to assess the amount of trapping from different mechanisms for long-term CO ₂ storage.

@en

We present an efficient compositional framework for simulation of CO2 storage in saline aquifers with complex geological geometries during a lifelong injection and migration process. To improve the computation efficiency, the general framework considers the essential hydrodynamic physics, including hysteresis, dissolution and capillarity, by means of parameterized space. The parameterization method translates physical models into parameterized spaces during an offline stage before simulation starts. Among them, the hysteresis behavior of constitutive relations is captured by the surfaces created from bounding and scanning curves, on which relative permeability and capillarity pressure are determined directly with a pair of saturation and turning point values. The new development also allows for simulation of realistic reservoir models with complex geological features. The numerical framework is validated by comparing simulation results obtained from the Cartesian-box and the converted corner-point grids of the same geometry, and it is applied to a field-scale reservoir eventually. For the benchmark problem, the CO2 is injected into a layered formation. Key processes such as accumulation of CO2 under capillarity barriers, gas breakthrough and dissolution, are well captured and agree with the results reported in literature. The roles of various physical effects and their interactions in CO2 trapping are investigated in a realistic reservoir model using the corner-point grid. It is found that dissolution of CO2 in brine occurs when CO2 and brine are in contact. The effect of residual saturation and hysteresis behavior can be captured by the proposed scanning curve surface in a robust way. The existence of capillarity causes less sharp CO2-brine interfaces by enhancing the imbibition of the brine behind the CO2 plume, which also increases the residual trapping. Moreover, the time-dependent characteristics of the trapping amount reveals the different time scales on which various trapping mechanisms (dissolution and residual) operate and the interplay. The novelty of the development is that essential physics for CO2 trapping are considered by the means of parameterized space. As it is implemented on corner-point grid geometries, it casts a promising approach to predict the migration of CO2 plume, and to assess the amount of CO2 trapped by different trapping mechanisms in realistic field-scale reservoirs.

@en

Various computational fluid dynamic simulations in engineering, such as external aerodynamics, only need the silhouette of an input geometry. Often, it is a laborious process that can take up many human hours. In addition, the CAD geometries are too complex and contain intricate features and topological holes. We showcase an effortless way to shrink-wrap triangulated surfaces with the sole intent of topology and surface simplification. Building upon the concepts of mathematical morphology and newer advancements in geometry processing, we present a straightforward and robust algorithm that can guarantee genus-zero surfaces. Our techniques are equally applicable to general polyhedral meshes and well-suited for handling both oriented and unoriented point clouds. We provide examples using unoriented point clouds to demonstrate the versatility of our algorithms. We have designed our algorithms with a wide variety of applications in mind. However, we specifically highlight their capability for aerodynamic simulations, fluid volume extraction, and surface simplification. Additionally, we emphasize the practicality and ease of implementing the proposed algorithms, and we chain additional algorithms to develop variants of our wrap algorithm.

@en

Constructing an analysis-suitable parameterization for the computational domain from its boundary representation plays a crucial role in the isogeometric design-through-analysis pipeline. PDE-based elliptic grid generation is an effective method for generating high-quality parameterizations with rapid convergence properties for the planar case. However, it may generate non-uniform grid lines, especially near the concave/convex parts of the boundary. In the present work, we introduce a novel scaled discretization of harmonic mappings in the Sobolev space H¹ to remit it. Analytical Jacobian matrices for the involved nonlinear equations are derived to accelerate the computation. To enhance the numerical stability and the speed of convergence, we propose a simple and yet effective preconditioned Anderson acceleration framework instead of using computationally expensive Newton-type iteration. Three preconditioning strategies are suggested, namely diagonal Jacobian, block-diagonal Jacobian, and full Jacobian. Furthermore, we discuss a delayed update strategy of the preconditioner, i.e., the preconditioner is updated every few steps to reduce the computational cost per iteration. Numerical experiments demonstrate the effectiveness and efficiency of our improved parameterization approach and the computational efficiency of our preconditioned Anderson acceleration scheme.

@en

Surrogate models based on convolutional neural networks (CNNs) for computational fluid dynamics (CFD) simulations are investigated. In particular, the flow field inside two-dimensional channels with a sudden expansion and an obstacle is predicted using an image representation of the geometry as the input. Generative adversarial neural networks (GANs) have been shown to excel at such image-to-image translation tasks. This motivates the focus of this work on investigating the specific effect of adversarial training on model performance. Numerical results show that the overall accuracy of the GANs is generally lower compared to an identical generator model trained directly on the ground truth using an L1 data loss. On the other hand, GAN predictions are often visually more convincing and exhibit a lower continuity residual.@en

The use of sequential time integration schemes becomes more and more the bottleneck within large-scale computations due to a stagnation of processor’s clock speeds. In this study, we combine the parallel-in-time Multigrid Reduction in Time method with a p-multigrid method to obtain a scalable solver specifically designed for Isogeometric Analysis. Numerical results obtained for two- and three-dimensional benchmark problems show the overall scalability of the proposed method on modern computer architectures and a significant improvement in terms of CPU timings compared to the use of standard spatial solvers.@en

In this article, we propose and analyze a set of fully nonstationary Anderson acceleration (AA) algorithms with two window sizes and optimized damping. Although AA has been used for decades to speed up nonlinear solvers in many applications, most authors are simply using and analyzing the stationary version of AA (sAA) with fixed window size and a constant damping factor. The behavior and potential of the nonstationary version of AA methods remain an open question. Most efficient linear solvers however use composable algorithmic components. Similar ideas can be used for AA to solve nonlinear systems. Thus in the present work, to develop nonstationary AA algorithms, we first propose a systematic way to dynamically alternate the window size (Formula presented.) by the multiplicative composite combination, which means we apply sAA((Formula presented.)) in the outer loop and apply sAA((Formula presented.)) in the inner loop. By doing this, significant gains can be achieved. Second, to make AA to be a fully nonstationary algorithm, we need to combine these strategies with our recent work on the nonstationary AA algorithm with optimized damping (AAoptD), which is another important direction of producing nonstationary AA and nice performance gains have been observed. Moreover, we also investigate the rate of convergence of these nonstationary AA methods under suitable assumptions. Finally, our numerical results show that some of these proposed nonstationary AA algorithms converge faster than the stationary sAA method and they may significantly reduce the storage and time to find the solution in many cases.

@en

We present the projection-based embedded discrete fracture model (pEDFM) for hexahedral corner-point grid (CPG) geometries, for the simulation of hydrothermal processes in fractured porous media. Unlike the previously-developed pEDFM for structured box grids, our new development allows for the modeling of complex geometries defined with hexahedral CPG cells. It also advances the pEDFM method to include coupled flow and heat transfer systems. Mass and energy conservation equations are simulated in a fully-coupled manner using a fully-implicit (FIM) integration scheme. This allows for stable simulations, specially when large time steps are taken. Independent corner-point grids are imposed on the rock matrix and all fractures, with conductivities ranging from highly permeable to flow barriers. The connectivities between the non-neighboring grid cells are described such that a consistent discrete representation of the embedded fractures occurs within the corner-point grid geometry, specially as the quadrilateral interfaces are allowed to be fully flexible. Various numerical tests including geologically-relevant and real-field models, which are established in the literature, are conducted to demonstrate the applicability of the developed method. It is shown that pEDFM can accurately capture the physical influence of both highly conductive fractures and flow barriers on the flow and heat transfer fields in complex reservoir geometries. This development is promising for flow simulations of real-field geo-models, increasing the discretization flexibility and enhancing the computational performance for capturing explicit fractures accurately.@en