To integrate renewable energies with our current energy systems, we require interaction between gas and electrical networks. The coupling of networks results in a larger system of equations to be solved. Henceforth, scalable solvers are more suitable for large coupled networks. In this paper, a preliminary research is done, by investigating Krylov solvers on gas networks from the GasLib library. The networks are simulated with steady-state models. The models yield a nonlinear system, which is solved with the Newton-Raphson method. The corresponding Jacobian is non-symmetric, indefinite and sparse. We have considered the following Krylov solvers: GMRES, Bi-CGSTAB and IDR(s). We compare the performance with a direct solver, which is the LU factorisation. Our results show that basic Krylov solvers are ineffective in solving the networks, because most networks have a large condition number and an unfavourable distribution of the eigenvalues. Hence, we have explored several preconditioners, such as Jacobi, Gauss-Seidel and ILU methods. Only the ILU preconditioner with the use of the COLAMD reordering scheme leads to convergence of all networks. For this preconditioner, the fill ratio has to be taken large enough, otherwise the ILU factorisation breaks down due to a zero pivot. The minimum required fill ratio leads to a similar amount of work as the direct solver. Thus the combination of ILU and Krylov solver does not perform better than direct solvers for these medium sized problems.

Large-scale geological storages of hydrogen (H₂) and carbon dioxide (CO₂) in saline aquifers present feasible options for a sustainable energy future. We compared the plume migration of CO₂ and H₂ in aquifers using the FluidFlower benchmark, incorporating the state-of-the-art thermophysical and petrophysical properties. The H₂ plume, with its higher buoyancy and mobility compared to CO₂, remains predominantly in the gas phase due to its lower solubility, increasing the chances of escaping through fractures or migration to distant regions. This additionally leads to a higher pressurized reservoir, which, along with higher buoyancy, increases the chance of caprock penetration. Dissolution trapping of CO₂ into brine increases over time due to its fingering, while H₂ does not show fingering. Our findings show that while geological carbon storage (GCS) benefits significantly from all structural, dissolution, and residual trapping, underground hydrogen storage (UHS) relies mainly on structural trapping, making the integrity of sealing elements of the system a key factor in its performance.

Infiltration models are crucial components of rainfall-runoff models based on shallow water equations, combined with direct-rainfall for flooding simulations. While the original Horton’s model is frequently used, various modifications have been proposed to deal with its limitations for intermittent rainfall patterns. We evaluate two modifications, from well-known Storm Water Management Model or SWMM and from Diskin and Nazimov, against the original Horton’s model both theoretically and with numerical experiments. We find the formulation from Diskin and Nazimov as most suitable for real-world applications. In this paper, we describe an adaptation of the Diskin and Nazimov infiltration model through a Surface Detention Box model to integrate it into a shallow water equation solver, that accounts for detailed topographic information using a sub-grid approach. The Surface Detention Box model in this paper is generalized to account for sources and sinks other than infiltration. We verify the efficiency of our implementation for a catchment in Australia with intermittent and extreme rainstorms. We also demonstrate the accuracy, efficiency and the precise volume conservation of our method for high-resolution grids and large computational time steps, enabled by the predictor–corrector solver. In conclusion, we present a robust and efficient scheme for practical flood simulations, including various sources and sinks such as rainfall and infiltration. Our approach is a strong foundation for operational flood forecasting with high resolution Digital Terrain Models.

We propose a linear electrolyzer model for steady-state load flow analysis of multi-carrier energy networks, where the electrolyzer is capable of producing hydrogen gas and heat. For our electrolyzer model, we show that there are boundary conditions that lead to a well-posed problem. We derive these conditions for two cases, namely with a known and unknown heat efficiency parameter. Furthermore, the derived conditions are validated numerically. Moreover, we investigate the extensibility of our model by including nonlinear models from electricity, gas, and heat. In this setting, we derived boundary conditions based on our previous findings. Due to the involvement of nonlinearity, it is a challenge to prove that the boundary conditions lead to a well-posed problem. Therefore, we simulated the electrolyzer connected with an electricity, gas, and heat system. Additionally, we considered a known and unknown heat efficiency parameter. The numerical results support that the linear electrolyzer model is solvable in a multi-carrier energy network.

The increasing reliance on 3D scanning and meshless methods highlights the need for algorithms optimized for point-cloud geometry representations in CAE simulations. While voxel-based binning methods are simple, they often compromise geometry and topology, particularly with coarse voxelizations. We propose an algorithm based on a Series of Local Triangulations (SOLT) as an intermediate representation for point clouds, enabling efficient upsampling and downsampling. This robust and straightforward approach preserves the integrity of point clouds, ensuring resampling without feature loss or topological distortions. The proposed techniques integrate seamlessly into existing engineering workflows, avoiding complex optimization or machine learning methods while delivering reliable, high-quality results for a large number of examples. Resampled point clouds produced by our method can be directly used for solving PDEs or as input for surface reconstruction algorithms. We demonstrate the effectiveness of this approach with examples from mechanically sampled point clouds and real-world 3D scans.

We present a matrix-free parallel scalable multilevel deflation preconditioned method for heterogeneous time-harmonic wave problems. Building on the higher-order deflation preconditioning proposed by Dwarka and Vuik (SIAM J. Sci. Comput. 42(2):A901-A928, 2020; J. Comput. Phys. 469:111327, 2022) for highly indefinite time-harmonic waves, we adapt these techniques for parallel implementation in the context of solving large-scale heterogeneous problems with minimal pollution error. Our proposed method integrates the Complex Shifted Laplacian preconditioner with deflation approaches. We employ higher-order deflation vectors and re-discretization schemes derived from the Galerkin coarsening approach for a matrix-free parallel implementation. We suggest a robust and efficient configuration of the matrix-free multilevel deflation method, which yields a close to wavenumber-independent convergence and good time efficiency. Numerical experiments demonstrate the effectiveness of our approach for increasingly complex model problems. The matrix-free implementation of the preconditioned Krylov subspace methods reduces memory consumption, and the parallel framework exhibits satisfactory parallel performance and weak parallel scalability. This work represents a significant step towards developing efficient, scalable, and parallel multilevel deflation preconditioning methods for large-scale real-world applications in wave propagation.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.