Research

We develop, distribute, and maintain peer-reviewed software with three purposes:

1. Disseminate our research on solvers and advanced processes in a reproducible manner.
2. Enable our industry partners to freely evaluate, utilize, reproduce, extend, and distribute, commercially or otherwise, our innovations.
3. Further our research into the science of how simulation software is developed, maintained, and distributed.

Our Software Research Mission

A quintessential target area for innovation today is to bridge the gap between mathematical semantics, physical notions, and computer programming languages. Owing to the physical complexity inherent to reservoir engineering problems, simulator developers have long recognized the pains and needs in this area. The goal of this research direction is to rethink and reengineer the way that we develop, deliver, and maintain simulation software in order to realize a future where:

1. Customization and ubiquitous interoperability are instantaneous and cost-free.
2. Domain-expert users are able to make rapid modifications to the Physics directly.

We aim to accomplish our mission as the culmination of three major directions:

• Discrete Problem Programming Toolkit DPPT
• Generic Solvers GenSOL
• Automated Simulator Generator System ASGS

Discrete Problem Programming Toolkit (DPPT)

A Discrete Problem is a software abstraction that implements the fully discretized equations that are solved by a simulator in order to obtain discrete approximations to the state variables in space and time.

As illustrated in the figures above, a Discrete Problem must offer a two-way interface that is well-engineered to admit a wide and growing plethora of numerical solution methods while simultaneously preventing the exposure of its internal particulars. A standardized interface will ensure effortless interchangeability of Discrete Problem implementations, thereby supporting both the Tailor- and General-Purpose paradigms.

At a minimum, a conforming Discrete Problem interface should be able to produce a discrete residual system and Jacobian matrix that are evaluated at given input states. Generic designs allow the parameterization of the Physics to any conforming numerical datastructures. The interface must also accommodate requests to adapt the internal specifics of the discrete numerical approximation.

Our research in this area focuses on developing a universal interface, and more importantly, on providing a software toolkit that cuts the development cost of any Discrete Problem to a minimum. This is to be accomplished by the automation of the activities entailed in developing a Discrete Problem.

The Discrete Problem Programming Toolkit (DPPT) is a multilayered library of highly configurable generic datastructures and algorithms. Modules within different layers can be used to develop a generic Discrete Problem using different semantics and levels of abstraction. Below is a depiction of the DPPT modules that are arranged in three layers.

1. Base Layer
The DPPT base layer offers several modules that implement efficient numerical data structures and algorithms. The semantics are closely tied to standard library concepts when efficiency is not sacrificed. For example, while standard library memory allocators cannot maintain state, our alternatives bypass these concept requirements in order to admit multithreaded pools. A beta point of departure for these modules is available in the faSTL library.

• Numerical datastructures
  • Dense collections including numeric arrays and vectors.
  • Indirection devices (iterators) that separate the storage pattern of variables in memory from their logical ordering.
  • Sparse collections including block sparse matrices and vectors.
  • Algebraic expression encodings using multiple implementations including expression templates.
• Memory managements
  • Heap allocators offering configurable pools and multithreading.
  • Casting and construction utilities.
• Thread Management
  • Interchangeable Policies for threading using multiple protocols including OpenMP, POSIX, and standard library protocols.
• Numerical algorithms
  • Polyalgorithmic sparse and dense linear combination kernels including recent algorithmic developments and expression templates.
  • Forced vectorization using inline SIMD instructions including SSE and intrinsics (platform dependent).
  • Miscellaneous: bit hacks, unrolling, and tiling .

2. Middle Layer
The middle layer offers modules that generically implement all aspects of discretizing and linearizing vector calculus operators. The semantics are mathematical. Configuration choices for all orthogonal aspects can be made by changing single tags within user code. For example, the following listing calculates the discrete nonlinear Laplacian and its analytical Jacobian matrix using single point upwinding and the MPFA O-method on a Cartesian topology that is ordered according to a Hilbert space-filling curve:

Variableset< memoryLayout, ADvector > U;DiscreteContinuum mesh( UniformCartesianBuilder< PeanoHilbert > build( “datafile.txt” ) );Laplacian< Omethod > ( mesh, U[TRANS].begin(), U[POTENTIAL].begin(), F.begin() );

• Discrete Operators: Implementations of finite element and finite volume approximations to nonlinear differential operators. These are independent of the mesh topology and variable layout in memory, and they are parameterized by the desired numerical scheme.
• Discrete Continuum: A generic representation of discrete spatial domains that may intersect. Alternate implementations are available in certain topologies for low-memory requirements in cases of massive meshes. These instances are built using input mesh data that are generated externally.
• Automatic Variables: Building on the initial success of the Automatically Differentiable Expression Templates Library (ADETL), we are continuing to develop a module that efficiently automates the analytical computation of derivatives of all expressions that occur in the callstack of a discrete Reservoir Simulation problem. The ADETL was designed to be extendable but it provides a limited scope of datastructures and algorithms to act as a proof of concept. Our work is to integrate and extend emerging AD algorithms and implementations to accommodate the wide range of expressions that are encountered in reservoir simulation.
  • Detailed profiling analysis and benchmarks for various stages in the computation of linearized equations.
  • Improved automatic activation of variables for variable switching and reordering.
  • Tailor approaches to univariate, multivariate dense, and multivariate sparse derivatives.
  • Facilities to convert types and to operate on arguments with heterogeneous expressions.
  • Polyalgorithmic forward and reverse mode AD with recommendations for appropriate use in different calculations.
  • Improved performance by anhialating expression construction costs.
  • Compressed AD using seed matrices to dramatically improve runtime costs in cases with predictable sparsity coutcomes.

3. Top Layer
The top layer module will offer abstractions for general partial differential equations such as the nonlinear transient advection-reaction-diffusion equations, Helmholtz and Poisson forms, and the time-dependent second-order wave equation.

GenSOL

We are designing and developing datastructure-neutral implementations of timestepping, nonlinear, and linear solution algorithms. The prime purpose for GenSOL is to provide a convenient platform with which to reproduce, analyze, compare, and even use all of our contributions in the solvers area.

Automated Simulator Generator System

Efforts such as our proposed DPPT promise to cut down on the costs of computer programming for the development of the physics portion of a simulator. Such efforts are exclusively tied to the underlying programming language. Moreover, they are typically only used by expert developers, barring the domain expert from direct participation. Problem Solving Environments such as “multiphysics packages” (e.g. COMSOL and the MRST) provide the means to develop physical models and simulation tools without the need for advanced programming language semantics. The efficiency and scalability of these tools however, remain far from ideal relative to compiled software. The purpose of the ASGS direction is simply to bridge this gap. We propose to develop a system that can translate domain-specific and high-level mathematical markup into C++ using the DPPT constructs. The schematic below illustrates this process.

The ASGS will ultimately enable cost-free customization of the physics within simulators. Domain experts can use a graphical equation editor to write, share, and disseminate “new models”, while vendors will be able to provide software that implements these models without incurring additional costs. The responsibility of devising well-posed physical models will be in the hands of the domain expert.

With growing physical complexity, the importance of nonlinear and linear solution processes is underscored. Our aim is to make fundamental advances in several challenge areas that arise ubiquitously in complex application scenarios. Leveraging these fundamental advances, we also aim to develop lasting, general, and robust solution technologies and methods. Our ultimate measure of success will be the wide and successful adoption of our methods in industrial and academic simulators to perform full-resolution, full-Physics modeling of seriously complex large-scale emerging processes.

Spatiotemporal Locality

Locality is inherent to all transient flow and transport phenomena. The superposition of the two disparate spatiotemporal characteristic scales that underly flow and transport leads to a problem of dimensionality. While the extent of an adequate simulation domain is dictated by the nature of flow, the resolution of its discrete form is restricted by the transport variables.

As a first response to this issue of dimensionality, Adaptive Discretization methods have been devised to exploit an a priori understanding of locality in order to reduce the computational cost of simulations. While such methods have provided various degrees of success in applications, it remains that they are inherently restricted by the fidelity of the discretization under aggressive adaptivity.

Our work in this areas seeks an alternate strategy. We are devising nonlinear solution processes that internally adapt the level of computation to precisely match that of the underlying spatiotemporal changes in state. This will enable the efficient use of high-fidelity, static, and full-resolution discretizations that do not require global computations unless absolutely necessary.

Fundamental Question: Over a Newton iteration, which cells will experience a nonzero update in state?

Encouraged by an initial result, we are answering this question for general flow and transport, independent of the spatial discretization that is used. The result is obtained by the application of the infinite dimensional Newton iteration process to the semidiscrete governing equations. That is, suppose that the governing equation has the form,

hen the semidiscrete fully-implicit form is simply,

In reservoir simulation, the generally nonlinear and heterogeneous flow term is discretized in space using an accurate method. We will not perform this discretization. rather we will work with the infinite dimensional (continuous) problem itself.

Over a single time-step, the initial state is,

and the subsequent state at the end of the time-step is,

Newton’s method can be applied to any differentiable nonlinear operator, including an infinite dimensional form. The iteration produces the sequence of iterations,

Note that this is a sequence of functions defined over the entire domain. The equation for the updates can be written as,

where the generalized Frechet derivative of the infinite dimensional residual equation is a linear partial differential equation. The Newton updates in a reservoir simulator are numerical approximations to the infinite dimensional update that can be obtained analytically in certain cases.

We have derived sharp conservative estimates for the analytical infinite dimensional Newton updates for general sequential flow and transport!

Using these estimates which are essentially free to evaluate, we can tell which cells within the domain will experience change over a Newton step, and which ones will not. A whitepaper describing this process is available here.

Our current research in this area is focused on:

1. Deriving estimates for the coupled form of reactive flow and transport in multiple dimensions.
2. Extensions to general nonlinear multicomponent advection, diffusion, and reaction problems including thermal.
3. Exploration of applications in fractured and stimulated media.

Application 1: solvers that only solve for cells that will experience change.

In any iterative nonlinear solution process that generates iterates that are the solution of a large sparse linear system, why solve for the entire system if we know that only a portion of the resulting vector has nonzero entries? Moreover, at initial iterates, far from the solution, it may be unnecessary to compute the linear solution very accurately. Our universal results on locality, provide an a priori identification of these vector components that is available at negligible cost. We are developing application grade solvers that do this. Compared to detection methods that are designed at the discrete level such as in this article, our methods are general to both flow and transport, and are available by evaluating a single formula!

Our current research in this area is focused on:

1. Devising solvers that do not degrade the convergence rate of the nonlinear iteration.
2. Developing strategies that adaptively alter the target linear convergence tolerance.

Application 2: Adaptive Mesh Refinement inside the nonlinear solver and not within the discretization.

Particularly in the case of fast physics such as flow, the nonlinear iteration updates may have relatively global support. In cases where these updates decay rapidly around certain points and very slowly near the boundaries, it is desirable to use adaptively localized meshes. Keeping the discretization mesh static, we are devising methods that solve the nonlinear residual using adaptive solver meshes. These ideas are inspired by nonlinear control optimization methods that adapt the resolution of control variables to reduce computational costs.

Application 3: Multirate nonlinear methods

Newton’s method is an explicit Euler discretization of the Newton Flow system of Ordinary Differential Equations using a stepsize of one. We propose to develop multirate or local timestepping methods to integrate the Newton Flow system. Using our results on spatiotemporal locality, we can decompose the spatial domain into several sets based on the anticipated rate of change in state.

Asymptotic Nonlinear Convergence

Safeguarding tactics are readily used to make Newton’s Method robust. That is, a safeguarded method is one that will ultimately converge to the answer if only it is given a chance to do so. This research direction builds upon recent advances in safeguarding to address the next question in the evolution of our understanding of converge.

Fundamental Question: How long will it take for a (safeguarded) Newton iteration to converge?

Our study on spatiotemporal locality is revealing. Newton updates for flow decay rapidly in the vicinity of nonzero entries in the residual. Away from such points, the decay is slow. The precise character of this decay depends on compressibility, time, and mobility. Newton updates for transport behave differently. The decay depends on an exponential ratio that is a nonlinear function of wave speeds.

Over one time-step, the number of iterations required for convergence is at least equal to the number of iterations required for the spatial support of non-zeros entries in the updates to encompass the number of cells that change over the time-step. We are deriving theoretical characterizations of this limit in terms of the decay parameters that have been derived to characterize locality.

This research will ultimately lead to an analogous result to the Courant-Friedrichs-Lewy condition for explicit time-step stability that is applicable to assessing the convergence rate of an implicit time-step using a safeguarded Newton iteration.

Since our results on locality are independent of the spatial discretization, it is conjectured that the convergence rate fo Newton is also independent of discretization. Rather, the important factors relate solely to the physics.

Application 1: Is it possible to somehow accelerate the nonlinear convergence rate of a safeguarded method?

Research is underway to develop nonlinear iterations that have an expanded spatiotemporal support. Our initial point of departure is to seek methods that are motivated by higher-order ODE integration. While expanded support has been observed with such methods, the additional computational costs of each step render them value-neutral. Current research is focused on retaining the expanded support properties while reducing the computational cost per iteration. Our approach to achieving this uses a differential geometry perspective, where the dynamics of motion can be abstracted with linear transformation operators. A promising idea is to develop extrapolation approaches on the Frenet Frame of the solution path.

Application 2: If we have bounds on the asymptotic convergence rates, we may better control timestep size selection for convergence.

A direct application of our results that characterize convergence rate is to use them in order to select timestep size in a manner that minimizes the computational cost per unit time advancement in implicit simulation. Suppose that at a specific timestep, time truncation error control dictates a maximum allowable timestep size,

Then, using our estimate for the number of iterations required for convergence as function of timestep size,

the objective is to select that target timestep size that satisfies the minimization problem,

Implicit Solvers With Embedded Timestep and Error Control

To date, reservoir simulation timestep selection strategies for error control are detached from the nonlinear solution process. Namely, there are two alternate approaches to timestep size selection for error control; a priori methods are exclusively based on extrapolation. Given state information from previous timesteps, estimates are extrapolated to constrain the next upcoming timestep size. On the other hand, a posteriori strategies will solve for a trial timestep and iterate to correct the trial timestep size until a more suitable one is attained.

Fundamental Question: Can we design implicit iterative solution methods that internally provide indications regarding time discretization errors?

Recently, a numerical Homotopty-Continuation procedure was proposed as an alternative to using Newton’s method to solve implicit timesteps. In this procedure, the homotopy parameter is the solution path arc-length. Iterations in this process solve jointly for a state and a corresponding timestep size. Subsequently, and unlike Newton’s method, each iterate is informative about the rate of change of the underlying physics in timestep. The mathematics of this alternate method are simple. At each implicit timestep, a nonlinear residual system of discrete equations is to be solved,

he method posits that for any positive timestep size, the corresponding state solution and the timestep size when augmented, form a point on a solution path. Because the solution path is a line, it may be parameterized by a single parameter such as arc-length or timestep itself,

and the initial point on the curve is simply the old state at the beginning of the timestep with corresponding timestep of zero,

Using the chain rule to differentiate the residual system with respect to the univariate parameter yields an Ordinary Differential Equation (ODE) whose solution is the solution path itself,

If the homotopy parameter is chosen as the timestep size itself, no additional equations are required for closure. More generally however, an additional equation is required, ensuring that,

By numerically integrating the ODE, the solution to any timestep size may be attained. Since each successive iteration is a solution to a successively increasing timestep size, we have the hope of embedding error control directly into this process. Start from a timestep size of zero, and proceed to integrate the ODE until local time discretization errors do not exceed the given desired threshold.

Ongoing Research: Numerical integration of the ODE system
Several facets of the required numerical integration procedure are currently under study. The overriding constraint to the design of a successful continuation method for implicit timestepping is that of computational efficiency.

As with any numerical integration process, there is drift; as the numerical integration proceeds, numerical approximation errors accumulate. Subsequently, continuation strategies are typically developed using either predictor-corrector approaches and/or adaptive discretizations. We are exploring the relative merits of such alternatives. To date, we have shown that it is possible to devise strategies with computational work demands that are on par with traditional Newton methods. Continuation methods that exceed the efficiency of safeguarded Newton methods on a timestep size basis have yet to be attained.

On the other hand, when the error control is used and/or the computational efficiency over the entire duration of a simulation are studied, the efficiency of continuation methods dominates.

Nonlinear Coupling of Multiphysics

There is growing interest in the large scale simulation of coupled multiphysics. The fully-coupled approach to solving such problems leads to a coupled nonlinear system. Our current understanding of nonlinear solution methods leads to the thought that different governing forms require different asymptotic convergence behaviors. For instance, while elliptic forms are unrestricted by the propagation dynamics, owing to the global spatial support of the linearized updates, parabolic and hyperbolic forms are much more so. Our research in this direction is to devise nonlinear solution methods that exploit this knowledge. Ultimately, components that are themselves invariant to propagation speeds can be updated less frequently than components that are wave-bound. The precise relative frequency of updates will depend on the level of coupling between the sets of physics. Our approach will focus on the infinite-dimensional forms.

Our ultimate goal is to leverage our work on flexible software and high-fidelity, efficient solvers towards the invention of future applications that are environmentally sound and economical.

A first-order target in this area is to contribute towards the ongoing evolution of our fundamental understanding of constitute relations such as those related to sorption, Knudsen effects in nanoscale zones, continuum damage, and multiphase flow. An initial example of such efforts is available here; simulation of CO2 sequestration and enhanced gas recovery in shale.

Beyond the first-order target, we also seek to experiment with ideas for groundbreaking (literally) applications using computational experimentation. We hope to marry fundamental research through our collaborations with field-scale process simulation and design.

Geomechanics

Stimulation methods such as hydraulic fracturing yield stimulated zones that are limited to a locale surrounding the well. Interactions with extensive natural fracture systems can lead to successful production from otherwise impermeable reservoirs. The interplay between pore pressure and geomechanical stresses and displacements is important. The effects could lead to dramatic differences in production as well as to questions regarding induced seismicity and/or structural stability. In addition to the accurate and practical simulation of steady-state geomechanical effects, we are exploring ideas that require us to incorporate time-dependent geomechanics.

Example: Lithotripsy is a non-invasive medical procedure that is used to break up stones in the kidney, bladder, or ureter. This process uses shock waves that are focused onto the target zone, causing stress concentrations that are sufficient to cause localized continuum damage.  We are investigating the potential application of this process to stimulate intrawell zones in the subsurface. After accurately coupling the stress wave propagation dynamics to flow and transport, we aim to investigate the impacts of wave dynamics (attenuation, refraction, and reflection), continuum damage models, and damage effects on flow properties. We will explore operating parameters including source locations and amplitudes in order to assess the feasibility of this process.

Electro-osmosis

Electro-osmosis is governed by a Helmholtz form. Under a voltage potential difference it is well-known that phases with different ionic polarity migrate to one pole or the other. Electro-osmosis is successfully applied in contaminant remediation. With recent interest in stimulation and recovery from shales, there are a number of potential applications in the energy sector.

Example: An important factor in drilling in shales is water migration from water based muds into the formation. The migration occurs due to a number of well-studied effects. Wellbore stability is frequently cited as a major concern that is compromised by water migration. A proposed idea is to apply electro-osmosis in order to reduce water migration out of the mud system and into the shale. By incorporating the governing form for electro-osmosis and coupling it to sorptive, advective, and diffusive transport, we hope to study the operating parameters that are required to improve drilling with water based muds in shales.

Coupled Free Flow and Transport

Processes such as Underground Coal Gasification involve the coupling between free mass and heat transport and flow in porous media. We aim to develop fully coupled models using well-established interface conditions at the boundaries.