The project team Pomdapi is concerned with the construction and analysis of simulation tools for the modeling of environmental and energy problems. These tools include numerical approximation schemes for partial differential equations, nonlinear solvers, and numerical techniques in optimization and complementarity problems. We are equally interested in reliable and correct programming methods for the implementation of these tools.
Most of the situations of interest in subsurface simulation involve more than one physical model. In these situations, one will need software modules for the simulation of the individual physical models, together with some model describing how these individual physics interact.
Geochemical phenomena result from the interaction of fluid and solute transport with chemical reactions in the earth's subsurface. These phenomena are important when modeling, for instance, a nuclear deep storage site, sequestration of carbon dioxide (CO2) in geological formations, or bioremediation of polluted sites.
The mathematical model for the coupled problem involves a multiphase flow model, describing the motion of several components and several phases through a porous medium, together with equations describing the various chemical reactions among these components. Coupling occurs through reactive source terms in the mass balance equations for the various components.
We have developed a coupling algorithm based on the Newton-Krylov method. It combines the strenghts of the more traditional operator splitting methods (separation of chemical and transport modules) with those of the more tightly coupled Newton's method (robustness, fast convergence).
The development of correct and efficient programs for scientific computing is a major challenge for the next decade.
We address the efficiency aspect by developping new algorithms and using parallelism techniques. We handle the correctness dimension at the programming level with the use of fixed set of templates or business objects, and at the theoretical level with the exploration of formal proofs of numerical programs.
It has been amply demonstrated that domain decomposition in space gives a scalable algorithm for elliptic (stationary) problems such as incompressible Darcy flow. The method can also be used for transport problems. In some of our applications (specifically in nuclear waste storage simulations), the physical properties in different regions of the domain vary over several orders of magnitude, making it desirable to use different discretization parameters (mesh size, and also time steps). Space-time domain decomposition methods make it possible to solve a transport problem using different time steps for different subdomains with different physics characteristics, such as permeabilities. The methods can even be formulated in the domain decomposition framework with non-matching grids in time and eventually also in space.
The most straightforward way to parallelize a code is to use a parallel library of solvers such as PETSc or Trilinos for iterative solvers, or MUMPS for direct solvers. These solvers are extremely efficient, and scalable. However, the black box approach may not always be suitable. As discussed in the previous section, domain decomposition methods have the potential of providing scalable solvers for a wide range of problems. They have enjoyed considerable success in a wide range of domains, ranging from solid mechanics to CFD. They can be implemented in a rather straightforward way on distributed memory computer, and are mostly insensitive to the type of grid used (though non-matching grids obviously require special care).
A common feature of both the physics and the computer architecture is their hierarchical nature: a coupled model (such as the reactive transport model seen above) involves several sets of equations, then a domain decomposition technique for solving these equations will introduce one additional level of parallelism. In most cases of interest the resulting subdomain problem may still be large enough to warrant a parallel solution, say by a parallel direct method.
The tools currently available for parallelization, while providing maximum flexibility, make it difficult to develop correct parallel programs. In contrast, the functional parallel skeleton programming approach offers a limited set of higher level constructs to define the parallel behavior of programs. We are implementing Sklml, a functional parallel skeleton compiler and programming system.
In the field of scientific computing, there is an increasing need for the implementation of algorithms that are more and more complex. For instance, the project-team deals with the following complex problems: (i) the coupling of numerical codes, including domain decomposition algorithms, and (ii) inverse problems, including optimal adaptive parameterization using refinement indicators. The applications solving those problems split into a generic algorithmic part, and a specific computing part. The idea is to provide once and for all the generic part as a driver executable program (the numerical business object); the user of the platform only needs to supply the worker programs specific to the application.
Our long-term goal of numerical program correctness is to obtain numerical codes that are entirely machine checked. This prospective research activity is a tight collaboration with the Epi Proval from Inria Saclay and LIPN from University of Paris 13, thanks to the ANR Fost. We already have formally proved a C program that implements a simple numerical scheme for the resolution of the one-dimensional acoustic wave equation. We now plan to tackle the formal proof of much more complex numerical programs including state-of-the-art programs using the finite element method.