Alternative dissolution-rate controlled model and time adaptive, high resolution scheme for site-scale subsurface carbon sequestration simulation

We present a parallelized geometric multigrid (GMG) method, based on the cell-based Vanka smoother, for higher order space-time finite element methods (STFEM) to the incompressible Navier–Stokes equations. The STFEM is implemented as a time marching scheme. The GMG solver is applied as a preconditioner for generalized minimal residual iterations. Its performance properties are demonstrated for 2D and 3D benchmarks of flow around a cylinder. The key ingredients of the GMG approach are the construction of the local Vanka smoother over all degrees of freedom in time of the respective subinterval and its efficient application. For this, data structures that store pre-computed cell inverses of the Jacobian for all hierarchical levels and require only a reasonable amount of memory overhead are generated. The GMG method is built for the deal.II finite element library. The concepts are flexible and can be transferred to similar software platforms.

During the past decade Graphics Processing Units (GPU) have been increasingly employed for speeding up compute intensive scientific applications. In this field, the geometric multigrid method (GMG) is one of the most efficient algorithms for solving large sparse linear systems of equations. Herein we analyze the performance of an optimized GPU based implementation of the GMG method on different state-of-the-art NVIDIA GPUs. The GTX Titan Black card, set-up with increased double precision performance leads to the smallest execution time. It is marginally faster than the more recently released GTX Titan X card which has considerably lower double precision performance. Moreover, an energy efficiency analysis reveals that the GTX 660M and the more powerful Titan cards require a similar amount of energy for running the GMG algorithm: the larger execution time is compensated by the lower power consumption.

A method is presented to include irregular domain boundaries in a geometric multigrid solver. Dirichlet boundary conditions can be imposed on an irregular boundary defined by a level set function. Our implementation employs quadtree/octree grids with adaptive refinement, a cell-centered discretization and pointwise smoothing. Boundary locations are determined at a subgrid resolution by performing line searches. For grid blocks near the interface, custom operator stencils are stored that take the interface into account. For grid block away from boundaries, a standard second-order accurate discretization is used. The convergence properties, robustness and computational cost of the method are illustrated with several test cases. New version program summaryProgram Title: AfivoCPC Library link to program files:https://doi.org/10.17632/5y43rjdmxd.2Developer's repository link:https://github.com/MD-CWI/afivoLicensing provisions: GPLv3Programming language: FortranJournal reference of previous version: Comput. Phys. Commun. 233 (2018) 156–166. https://doi.org/10.1016/j.cpc.2018.06.018Does the new version supersede the previous version?: Yes.Reasons for the new version: Add support for internal boundaries in the geometric multigrid solver.Summary of revisions: The geometric multigrid solver was generalized in several ways: a coarse grid solver from the Hypre library is used, operator stencils are now stored per grid block, and methods for including boundaries via a level set function were added.Nature of problem: The goal is to solve Poisson's equation in the presence of irregular boundaries that are not aligned with the computational grid. It is assumed these irregular boundaries are defined by a level set function, and that a Dirichlet type boundary condition is applied. The main applications are 2D and 3D simulations with octree-based adaptive mesh refinement, in which the mesh frequently changes but the irregular boundaries do not.Solution method: A geometric multigrid method compatible with octree grids is developed, using a cell-centered discretization and point-wise smoothing. Near irregular boundaries, custom operator stencils are stored. Line searches are performed to locate interfaces with sub-grid resolution. To increase the methods robustness, this line search is modified on coarse grids if boundaries are otherwise not resolved. The multigrid solver uses OpenMP parallelization.

The purpose of this paper is to discuss parallel preconditioning techniques to solve the elliptic portion (since it dominates computation) of the bidomain model, a non-linear system of partial differential equations that is widely used for describing electrical activity in the heart. Specifically, we assessed the performance of parallel multigrid preconditioners for a conjugate gradient solver. We compared two different approaches: the Geometric and Algebraic Multigrid Methods. The implementation is based on the PETSc library and we reported results for a 6-node Athlon 64 cluster. The results suggest that the algebraic multigrid preconditioner performs better than the geometric multigrid method for the cardiac bidomain equations.KeywordsCoarse GridMultigrid MethodParallel SpeedupConjugate Gradient IterationParallel GeometricThese keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

Although Gaussian elimination with partial pivoting is a robust algorithm to solve unsymmetric sparse linear systems of equations, it is difficult to implement efficiently on parallel machines because of its dynamic and somewhat unpredictable way of generating work and intermediate results at run time. In this paper, we present an efficient parallel algorithm that overcomes this difficulty. The high performance of our algorithm is achieved through (1) using a graph reduction technique and a supernode-panel computational kernel for high single processor utilization, and (2) scheduling two types of parallel tasks for a high level of concurrency. One such task is factoring the independent panels in the disjoint subtrees of the column elimination tree of A. Another task is updating a panel by previously computed supernodes. A scheduler assigns tasks to free processors dynamically and facilitates the smooth transition between the two types of parallel tasks. No global synchronization is used in the algorithm. The algorithm is well suited for shared memory machines (SMP) with a modest number of processors. We demonstrate 4- to 7-fold speedups on a range of 8 processor SMPs, and more on larger SMPs. One realistic problem arising from a 3-D flow calculation achieves factorization rates of 1.0, 2.5, 0.8, and 0.8 gigaflops on the 12 processor Power Challenge, 8 processor Cray C90, 16 processor Cray J90, and 8 processor AlphaServer 8400.

An oil/water, 3-D fully implicit Reservoir Simulator has been implemented on the ICL Distributed Array Processor (DAP), by a joint BNOC/ICL project team. The DAP is a single instruction, multiple data stream machine with 4096 processors in parallel. The major computational areas of code have been implemented on the DAP, leaving I/O and well calculations as serial code on the host 2976 machine. Matrix assembly is inherently parallel and can be implemented easily and efficiently on the DAP with compact code. Table look-up, though not obviously parallel, has been coded using a very efficient parallel algorithm. The Linear Solver, of Line Gauss-Seidel type, has been chosen for its efficiency on the DAP, with odd and even reduction used for solution of the resulting tri-diagonal equations. Further research will be undertaken to increase the robustness of this solver. The initial implementation has a 1-1 mapping of active grid blocks to processors, thus allowing up to 4096 active grid blocks. It is planned to extend to much larger models, and to 3 phase.

In this paper three different schemes for the finite-difference approximation of ( ) terms with irregular spacings of the grid points are investigated. Several methods of selecting the grid spacing are considered, and numerical results are presented for some carefully selected examples. presented for some carefully selected examples. Results indicate that one of the most commonly used finite-difference schemes may unnecessarily result in large discretization errors; the other two schemes discussed result in greater accuracy and are no more difficult to use. Introduction Finite-difference methods are used widely for solving partial-differential equations encountered in petroleum reservoir simulation. In many practical problems it is necessary to refine the grid in certain problems it is necessary to refine the grid in certain parts of the reservoir, thus resulting in an irregular parts of the reservoir, thus resulting in an irregular grid spacing for the system. For example, the local refinement of grid is usually necessary in the study of the flow behavior near a well, i.e., the case of coning problems. On the other hand, it is usually advantageous to make the grid locally coarse over aquifers and large gas caps. In spite of the importance of irregular grids in reservoir simulation and other engineering problems, little is known about the effect of various schemes on the over-all discretization error. Most literature on the subject of irregular grid spacings is concerned with the treatment of irregular boundaries. Here, the properties of several difference approximations to properties of several difference approximations to the second-order operator Au = [ ] are investigated. Operator A is encountered in all diffusion-type equations, such as multiphase flow in porous media and heat conduction with temperature-dependent conductivity. For simplicity, we limit our study to one space variable in this study. All results reported here can be extended readily to any number of space variables.

In this article we present a complete parallelization approach for simulations of PDEs with applications in optimization and uncertainty quantification. The method of choice for linear or nonlinear elliptic or parabolic problems is the geometric multigrid method since it can achieve optimal (linear) complexity in terms of degrees of freedom, and it can be combined with adaptive refinement strategies in order to find the minimal number of degrees of freedom. This optimal solver is parallelized such that weak and strong scaling is possible for extreme scale HPC architectures. For the space parallelization of the multigrid method we use a tree based approach that allows for an adaptive grid refinement and online load balancing. Parallelization in time is achieved by SDC/ISDC or a space-time formulation. As an example we consider the permeation through human skin which serves as a diffusion model problem where aspects of shape optimization, uncertainty quantification as well as sensitivity to geometry and material parameters are studied. All methods are developed and tested in the UG4 library.

The Geometric Multigrid (GMG) method is widely used in numerical analysis to accelerate the convergence of partial differential equations solvers using a hierarchy of grid discretizations. Multiple grid sizes and recursive expression of multigrid cycles make the task of program optimization tedious. A high-level language that aids domain experts for GMG with effective optimization and parallelization support is thus valuable.

Sparse approximate inverse smoothers for geometric and algebraic multigrid

Novel high resolution schemes for convective flow approximation are developed and coupled with general continuous full-tensor Darcy flux approximations. This development leads to new locally conservative formulations for multiphase flow in porous media. The higher order schemes are designed for fluid transport on unstructured grids and are constructed using slope limiters such that they are stable with a maximum principle that ensures solutions are free of spurious oscillations. The schemes are developed to handle unstructured meshes with variable grid spacing and different formulations are compared. Benefits of the resulting schemes are demonstrated for classical test problems in reservoir simulation including cases with full tensor permeability fields. The test cases involve a range of unstructured grids with variations in grid spacing, orientation and permeability that lead to flow fields that are poorly resolved by standard simulation methods. The new formulation is compared with standard reservoir simulation schemes and control-volume finite element methods. The new schemes are shown to effectively reduce numerical diffusion while resolving flow induced by rapid changes in permeability, leading to superior resolution of concentration and saturation fronts compared to other schemes.

Local uniform grid refinement (LUGR) is an adaptive grid technique for computing solutions of partial differential equations possessing sharp spatial transitions. Using nested, finer-and-finer uniform subgrids, the LUGR technique refines the space grid locally around these transitions, so as to avoid discretization on a very fine grid covering the entire physical domain. This paper examines the LUGR technique for time-dependent problems when combined with static regridding. Static regridding means that in the course of the time evolution, the space grid is adapted at discrete times. The present paper considers the general class of Runge-Kutta methods for the numerical time integration. Following the method of lines approach, we develop a mathematical framework for the general Runge-Kutta LUGR method applied to multispace-dimensional problems. We hereby focus on parabolic problems, but a considerable part of the examination applies to hyperbolic problems as well. Much attention is paid to the local error analysis. The central issue here is a "refinement condition" which is to underly the refinement strategy. By obeying this condition, spatial interpolation errors are controlled in a manner that the spatial accuracy obtained is comparable to the spatial accuracy on the finest grid if this grid would be used without any adaptation. A diagonally implicit Runge-Kutta method is discussed for illustration purposes, both theoretically and numerically.

A Survey of the Parallelization Issues for Geometric and Algebraic Multigrid Methods based on Generic Banded Approximate Inverses

In 1945 Hassler described a method for relative permeability measurement that was conceptually permeability measurement that was conceptually well founded but from the start has presented laboratory workers with various measurement problems. Recent as well as earlier analyses of problems. Recent as well as earlier analyses of the situation, for example, have been motivated by the need for an independent standard against which relative permeability data obtained by other (more convenient, less exacting) techniques can be calibrated. The latter is wanted as a hedge against the possibility that reservoir simulation predictions of possibility that reservoir simulation predictions of reserves and recovery will exhibit marked sensitivity to the input relative permeability data.The purpose of this paper is to call attention to a previously unmentioned type of measurement previously unmentioned type of measurement problem that appears to be inherent to the Hassler problem that appears to be inherent to the Hassler flow systems (Fig. 1). For simplicity, attention will be limited to two-phase colinear flow of a wetting and a nonwetting fluid. In such a case, the intention is that the pressure drop be identical in the contiguous immiscible fluids so that the capillary pressure gradient (and, hence, the saturation gradient) will be zero along the flow path when the relative permeability measurements are being made. permeability measurements are being made.For example, suppose that in a particular run a pressure, p1, is set where the wetting fluid enters the pressure, p1, is set where the wetting fluid enters the inflow capillary barrier, and a pressure, p4 less than p1, is set where the wetting fluid exits from the outflow capillary barrier. Then, as schematized in Fig. 2, the internal wetting-fluid pressures at the surfaces of discontinuity between the core sample and the capillary barriers, p2 and p3 less than p2, will assume values according to Darcy's law that depend on factors such as the ratio of the length of the core sample to that of an end barrier M and the ratio of the effective permeability of the core (to the flow of the wetting permeability of the core (to the flow of the wetting fluid) to the absolute permeability of the barrier material N. In fact, it is easy to show that the relationships under discussion are (1a) and (1b) The critically important things in Hassler work, of course, is to have the pressures p2 and p3 measured or otherwise immediately known just as soon as two-phase flow is being initiated at each step in the measurement of the relative permeability curves. This is because the relative permeability data are to be referred to particular values of saturation and types (e.g., drainage, imbibition) of saturation distribution and specifically to those that are prescribed by the particular values of capillary prescribed by the particular values of capillary pressure, Pc, that have been imposed before each pressure, Pc, that have been imposed before each two-phase flow episode. In other words, the inflow and outflow pressures in the nonwetting fluid, p5 and p6, respectively, must be set at precisely the same time p1 and p4 are set (i.e., as two-phase flow is initiated simultaneously in both fluids) in such a way that the following condition is met. (2) The dashed curves of Fig. 2, for example, depict how pressure would be distributed in the wetting fluid for a representative case where the constraint condition, Eq. 2, has not been met. The consequence of encountering such a situation in a Hassler-like relative permeability run is illustrated in Figs. 3 and 4. The dashed curves indicate that a positive capillary pressure gradient gives rise to a negative saturation pressure gradient gives rise to a negative saturation gradient in the direction of two-phase flow (and vice versa). P. 1161

Salinization of coastal aquifers is a current problem in many regions worldwide. Simulation of this process provides an important tool for the forecast of drinking water resources. The consideration of the unsaturated phreatic zone has an essential influence on the accuracy of the prediction. However, for the large temporal and spatial scales of the aquifers, the presence of the unsaturated subdomains creates challenging difficulties for the numerical methods. In this work, we investigate two approaches for simulating haline density-driven flow in partially saturated aquifers. The first approach, based on the Richards equation, uses a representation of the saturation field and employs an adaptive linearly implicit time discretization scheme. The second approach explicitly represents the water table using a level-set method and relies on weak coupling combined with a standard implicit Euler scheme. Both approaches employ a finite-volume discretization for spatial representation of fluid flow and salt transport. The implementation is based on the UG4 toolkit together with the parallel groundwater flow simulator d3f++. The resulting linear systems are solved using a geometric multigrid method. We compare the aforementioned approaches with respect to theoretical and numerical aspects. The performance of both methods is evaluated in three numerical experiments. In particular we demonstrate robustness of the linearly implicit scheme and introduce a stabilization for the level-set method. The high-performance computing (HPC) potential of the approaches is assessed and demonstrated as well.