Search

International audience; The Circulant Embedding Method (CEM) is a well known technique to generate stationary Gaussian Random Fields (GRF). The main idea is to embed the covariance matrix in a larger nested block circulant matrix, whose factorization can be rapidly computed thanks to the fast Fourier transform (FFT) algorithm. The CEM requires the extended matrix to be at least positive semidefinite which is proven to be the case if the enclosing domain is sufficiently large, as proven by Theorem 2.3 in [Graham et al, SIAM Journal on Numerical Analysis, 2018] for cubic domains. In this paper, we generalize this theorem to the case of rectangular parallelepipeds. Then we propose a new initialization stage of the CEM algorithm that makes it possible to quickly jump to a domain size close to the one needed for the CEM algorithm to work. These domain size estimates are based on fitting functions. Examples of fitting functions are given for the Mat\'ern family of covariances. These functions are inspired by our numerical simulations and by the theoretical work from [Graham et al, SIAM Journal on Numerical Analysis, 2018]. The parameters estimation of the fitting functions is done numerically. Several numerical tests are performed to show the efficiency of the proposed algorithms, for both isotropic and anisotropic Matérn covariances.

International audience; Reactive transport plays an important role in various subsurface applications, including carbon dioxide sequestration, nuclear waste storage, biogeochemistry and the simulation of hydro-thermal reservoirs. The model couples a set of partial differential equations, describing the transport of chemical species, to nonlinear algebraic or differential equations, describing the chemical reactions. Solution methods for the resulting large nonlinear system can be either fully coupled or can iterate between transport and chemistry. This paper extends previous work by the authors where an approach based on the Newton-Krylov method applied to a reduced system has been developed. The main feature of the approach is to solve the nonlinear system in a fully coupled manner while keeping transport and chemistry modules separate. Here we extend the method in two directions. First, we take into account mineral precipitation and dissolution reactions by using an interior point Newton method, so as to avoid the usual combinatorial approach. Second, we study two-dimensional heterogeneous geometries. We show how the method can make use of an existing transport solver, used as a black box. We detail the methods and algorithms for the individual modules, and for the coupling step. We show the performance of the method on synthetic examples.

International audience This paper develops a posteriori estimates for global-in-time, nonoverlapping domain decomposition (DD) methods for heterogeneous diffusion problems. The method uses optimized Schwarz waveform relaxation (OSWR) with Robin transmission conditions on the space-time interface between subdomains, and a lowest-order Raviart–Thomas–Nédélec discretization in the subdomains. Our estimates yield a guaranteed and fully computable upper bound on the error measured in the space-time energy norm of [19, 20], at each iteration of the space-time DD algorithm, where the spatial discretization, the time discretization, and the domain decomposition error components are estimated separately. Thus, an adaptive space-time DD algorithm is proposed, wherein the iterations are stopped when the domain decomposition error does not affect significantly the global error, allowing important savings in terms of the number of domain decomposition iterations while guaranteeing a user-given precision. Numerical results for a two-dimensional heat equation are presented to illustrate the efficiency of our a posteriori estimates and the performance of the adaptive stopping criteria for the space-time DD algorithm.

International audience; In this paper we are interested in the "fast path" fracture and we aim to use global-in-time, nonoverlapping domain decomposition methods to model flow and transport problems in a porous medium containing such a fracture. We consider a reduced model in which the fracture is treated as an interface between the two subdomains. Two domain decomposition methods are considered: one uses the time-dependent SteklovPoincaré operator and the other uses optimized Schwarz waveform relaxation (OSWR) based on Ventcell transmission conditions. For each method, a mixed formulation of an interface problem on the space-time interface is derived, and different time grids are employed to adapt to different time scales in the subdomains and in the fracture. Demonstrations of the well-posedness of the Ventcell subdomain problems is given for the mixed formulation. An analysis for the convergence factor of the OSWR algorithm is given in the case with fractures to compute the optimized parameters. Numerical results for two-dimensional problems with strong heterogeneities are presented to illustrate the performance of the two methods.

International audience; The far field simulation of underground nuclear waste disposal site requires a high computational cost due to the widely varying properties of different materials, the different length and time scales, and the high accuracy requirements. Nonoverlapping domain decomposition methods allow local adaptation in both space and time and result in parallel algorithms. We have extended the optimized Schwarz waveform relaxation (OSWR) method, successfully used for finite elements and finite volumes, to the case of mixed finite elements with their local mass-conservation property. Another choice is the substructuring method, which has been shown to be efficient for steady state problems with strong heterogeneities. We study a time-dependent Schur complement method, which is the algebraic counterpart of the discrete Steklov Poincaré operator, and introduce the Neumann preconditioner as well as weight matrices designed to make the convergence speed independent of the heterogeneities. Both methods enable the use of local time steps when the subdomains have highly different physical properties. Their performance is illustrated on test cases suggested by nuclear waste disposal problems.

National audience; Notre objectif est d'étudier la modélisation de phénomènes couplés et leur mise en oeuvre sur une grille de calcul, en développant quatre applications. Les deux premiers couplages sont de type physico-chimique et physico-physique (couplage algébrique d'équations). Ils sont appliqués respectivement à la contamination d'aquifères et à l'intrusion d'eau salée. Les deux autres couplages sont géométriques de type multidomaines ou multi-échelles et sont appliqués au stockage profond de déchets radioactifs, d'une part dans un milieu peu fracturé, d'autre part dans un réseau de fractures.Nous choisissons une approche par composants logiciels, qui permet d'encapsuler chaque code modélisant un phénomène physique. Les interfaces des composants permettent d'eectuer les échanges de données nécessaires au couplage numérique.L'exécutif PadicoTM garantit un calcul à haute performances sur une grille de calcul avec différents types de réseaux, grâce notamment à un modèle de composants parallèles, qui permet de passer à l'échelle.

Andra, CEA and EDF are conjointly developing a software platform for waste disposal simulations. The core of the tool for solute plume evaluation is the geochemistry transport coupling algorithm. Within Alliances, two geochemistry codes, PhreeqC and Chess, and three transport codes, Castem, Traces and Mt3D are available and interact within a sequential iterative algorithm for geochemistry transport coupling. Considering the waste disposal, one of the main problematic of a repository is to evaluate the alkaline perturbation in the host rock due to the cement interface with bentonite and clay. In order to deal with the heterogeneities and the various scales of the problem in two dimensions of space, a mesh adaptation method has been introduced to better deal with sharp interfaces. The adaptation algorithm will be described in details, pointing out its specificities due to the geochemistry treatment, especially concerning mineral interfaces. Numerical results will illustrate the algorithm performance for convection dominant and diffusion dominant problems.

In complex systems the domain of calculation for flow and contaminant transport may be the union of several subdomains with very different hydrogeological properties. Because of these differences, the subdomains have different space-time scales. One can think of adjusting the time scales with an almost uniform mesh to obtain a nearly uniform error, but an alternative is to use different space discretizations with a unique time discretization. This implies that the subdomains will have to be discretized independently and consequently the subdomain meshes will not match at the interfaces between subdomains. In order to solve this difficulty mortars and nonmortars techniques were developed. We use here a technique with mortars, that is a technique introducing extra unknowns on the interface between subdomains. This technique is combined with mixed finite elements and the resulting linear system is solved with a domain decomposition method along the lines of [1]. An example from the study of the flow and transport around a deep underground nuclear waste disposal will show the possibilities of the method. [1] G. Pencheva and Y. Yotov, Balancing decomposition for mortar mixed finite element methods, Numerical Algebra With Applications 10 (2003), pp. 159-180.

The simulation of multispecies reacting systems in porous media is of importance in several different fields: for computing the near field in nuclear waste simulations, in the treatment of bio-remediation, and in the evaluation of underground water quality. When looking at a coupled system, the role of chemistry is simply to separate the species into mobile and immobile species (immobile species may come from either sorption or precipitation reactions, and are not subject to transport). In this work, we assume that the medium is saturated, and that surface reactions do not change the porosity. Multi-species chemistry involves the solution of ordinary differential equations (if the reactions are kinetic) or nonlinear algebraic equations (if we assume local equilibrium). When simulating a coupled system, these equations have to be solved at each (grid) point, and at every time (step), leading to a huge coupled non-linear system. As has been observed several times, it is essential to use efficient numerical methods so as to be able to handle the size of systems occurring in the applications. If we assume local equilibrium, than the coupled system may be written as a DAE. Since methods and software for solving DAEs have reached a high level of maturity (at least for small index systems, which is the case here), it is natural to try and use this technology. One important issue is the size of the system to be solved, as all chemical species at all grid points are coupled. For any realistic configuration, it will not be possible to form, let alone factor, the Jacobian matrix. A better solution is to use Newton--Krylov methods, where the linear system at each Newton iteration is solved by an iterative method. We can thus keep the fast convergence of Newton's method, while only requiring Jacobian matrix--vector products, and these can be approximated by finite differences. We will discuss two methods, based on different choices for the main unknowns: In the ``one level'' method, we take as unknowns both the total dissolved and fixed concentrations, but also the concentrations for the component species that are the main unknowns of the chemical problem. This is closely related to the Direct Substitution Approach used by the geochemists. By carefully formulating the problem, we can put it into the DAE framework. At each time step, a single non-linear system has to be solved, and the Jacobian matrix couples both the transport matrix, and the chemistry Jacobian matrix. The method is called one-level as we solve transport and chemistry at the same time. Provided that we can access this Jacobian matrices (and this may not always be possible in practice), this may be less expensive than the two-level methods below. The second, ``two-level'' method takes as unknowns only the total dissolved and fixed concentrations. The system to be solved at each time step requires the solution of the chemical sub-problem at each Newton iteration, when evaluating the residual. On the orher hand, the Jacobian matrix again involves the transport matrix, as well as the Jacobain of the ``chemical solution operator'' mapping. This is only known implicitly, and the Jacobian has again to be evaluated by finite diffferences. Even though this method may be more expensive than the one-level method, its main advantage is to make it possible to treat chemistry as a black--box, even in the Newton--Krylov context. This may be important, as chemical simulators are becoming increasingly sophisticated. A major issue for both methods will be the choice of a suitable preconditioner.

In this presentation we are concerned with the problem of developing a numerical simulator for modeling flow and contaminant transport in the vicinity of an underground nuclear waste repository for use in performance assessment. Clearly reasonably accurate and reliable simulators are needed, but as performance assessment requires many simulations, it is also important that the simulators be very efficient. The simulations under consideration necessarily involve large scale calculations both in space and in time, and yet must take into account important heterogeneities in the domain. This suggests the use of multidomain modeling: the domain is divided into more nearly homogeneous subdomains and efficient numerical solvers are used for each subdomain. The calculations on the various subdomains are related through transmission or interface conditions, and parallel computing and high level programming languages can then be used for numerical implementation. Around the waste depository the porous medium is made up of several regions of different permeability, the permeability varying over several orders of magnitude from one region to the next. In the method developed here the domain is partitioned into homogeneous subdomains and nonoverlapping domain decomposition is used for both the flow calculations and to solve the transport equation. The single phase flow equation is an elliptic equation coupling a conservation equation with Darcy's law. As the Darcy velocity is the unknown of interest, a mixed finite element method is used to solve the equation within the subdomains. In order to be able to use general hexahedral grids, Kuznetsov-Repin macroelements are used. Continuity of the pressure and continuity of the flux are enforced at the interfaces between the subdomains using a balancing domain decomposition method. This method has the advantage of being very robust with respect to large contrasts in permeability. The transport equation is of advection-diffusion type. To approximate the solution to this equation an operator splitting technique is used: an advection equation is solved alternately with a diffusion-dispersion-decay equation. The parabolic diffusion-dispersion-decay equation is solved implicitly using mixed finite elements while the advection equation is solved explicitly using cell centered finite volumes with upstream weighting. In this way we are able to use a smaller time step for the advection equation when the advection coefficient is large. The domain decomposition methods used for the pressure equation extend also to these time dependent problems. In subdomains of high permeability, the natural time scale is much more rapid than in neighboring subdomains with lower permeability, and a much smaller time step size must be used for the advection step is these subdomains. However if only one time step size is used for all of the advection steps, it is the more rapid time scale that is imposed and the much smaller time step size must be used throughout the domain even if the size of the region in which it is required is relatively small. This suggests the use of different time step sizes for different subdomains. One of the many advantages of using domain decomposition techniques is that they lend themselves readily to parallelization. However, implementing the subdomain coupling requires complex programming. This is done efficiently using the programming environment OCamlP3l. This is a recently developed part of the strongly typed functional programming language OCaml which can be used for parallel computing. Realistic three dimensional numerical experiments from the Andra project will be shown.