David Moulton › Publications

Modeling of multiphase flow and transport in highly heterogeneous porous media must capture a broad range of coupled spatial and temporal scales. Recently, a hierarchical approach dubbed the Multilevel Multiscale Mimetic (M³) method, was developed to simulate two-phase flows in porous media. The M³ method is locally mass conserving at all levels in its hierarchy, it supports unstructured polyhedral grids and full tensor permeabilities, and it can achieve large coarsening factors. In this work we consider infiltration of water into a two-dimensional layered medium. The grid is aligned with the layers but not the coordinate axes. We demonstrate that with an efficient temporal updating strategy for the coarsening parameters, fine-scale accuracy of prominent features in the flow is maintained by the M³ method.

Black Box Multigrid (BoxMG) is a robust variational multigrid solver for diffusion equations on logically structured grids. BoxMG standardly uses coarsening by a factor of two. It handles cell-centered discretizations on logically rectangular grids by treating the cell-centers as the unknowns to be coarsened. Such a strategy does not preserve the cell structure. That is, coarse-grid cells are not the union of fine grid cells. In some applications, such as local grid refinement, it is desirable that the cell structure be preserved. In this paper we develop a method that employs coarsening by a factor of three. It is a natural generalization of standard BoxMG, using operator-induced interpolation (which approximately preserves the continuity of the normal flux), restriction as the transpose of interpolation, and Galerkin coarsening. We present numerical results that demonstrate its robustness with respect to discontinuous diffusion coefficients.

We describe a multilevel multiscale mimetic (M³) method for solving two-phase flow (water and oil) in a heterogeneous reservoir. The governing equations are the elliptic equation for the reservoir pressure and the hyperbolic equation for the water saturation. On each time step, we first solve the pressure equation and then use the computed flux in an explicit upwind finite volume method to update the saturation. To reduce the computational cost, the pressure equation is solved on a much coarser grid than the saturation equation. The coarse-grid pressure discretization captures the influence of multiple scales via the subgrid modeling technique for single-phase flow recently proposed in [?]. We extend significantly the applicability of this technique by developing a new robust and efficient method for estimating the flux coarsening parameters. Specifically, with this advance the M³ method can handle full permeability tensors and general coarsening strategies, which may generate polygonal meshes on the coarse grid. These problem dependent coarsening parameters also play a critical role in the interpolation of the flux, and hence, in the advection of saturation for two-phase flow. Numerical experiments for two-phase flow in highly heterogeneous permeability fields, including layer 68 of the SPE Tenth Comparative Solution Project, demonstrate that the M³ method retains good accuracy for high coarsening factors in both directions, up to 64 for the considered models. Moreover, we demonstrate that with a simple and efficient temporal updating strategy for the coarsening parameters, we achieve accuracy comparable to the fine-scale solution, but at a fraction of the cost.

The Laplace-Beltrami system of nonlinear, elliptic, partial differential equations has utility in the generation of computational grids on complex and highly curved geometry. Discretization of this system with Finite Elements readily accommodates unstructured grids, but generates a large, sparse, ill-conditioned system of nonlinear discrete equations. The extensive use of the Laplace-Beltrami approach, particularly in large-scale applications, has been limited by the scalability and efficiency of solvers. This paper addresses this limitation by developing two nonlinear solvers based on the Jacobian-Free Newton-Krylov (JFNK) methodology. A key feature of these methods is that the Jacobian is not formed explicitly for use by the underlying linear solver. Iterative linear solvers such as the Generalized Minimal RESidual (GMRES) method do not technically require the stand-alone Jacobian; instead its action on a vector is approximated through two nonlinear function evaluations. The preconditioning required by GMRES is also discussed; two different preconditioners are developed, both of which are readily treated with existing Algebraic Multigrid (AMG) methods. Further, the most efficient preconditioner overall for the problems considered is based on a Picard linearization. Numerical examples demonstrate that these new solvers are significantly faster than a standard Newton-Krylov approach; a speedup factor of approximately 26 was obtained for the Picard preconditioner on the largest grids studied here. In addition, these JFNK solvers exhibit good algorithmic scaling with increasing grid size.

In 2005, Christen and Fox introduced a delayed acceptance Metropolis-Hastings (DAMH) algorithm that improved computational efficiency in sample-based imaging of electrical conductivity (EIT). That work used a linear approximation to the forward map in the first step of the algorithm. In this paper, we develop an alternative approximation for use in DAMH, namely a multilevel approximation developed from the hierarchy of coarse-scale models obtained by variational coarsening. This approach builds on two important strengths of robust multigrid solvers. First, the cost of a fine-scale solution of the forward map scales linearly with the degrees of freedom, and hence, it is provides better efficiency for algorithms performing sample-based inference. Second, the homogenization implicit in robust variational multigrid methods gives better solutions at coarse scales than homogenization by averaging of coefficients. We report results from a stylized example in electrical impedance imaging where data is a noisy and incomplete measurement of the Dirichlet-to-Neumann map.

In 2005, Christen and Fox introduced a delayed acceptance Metropolis-Hastings (DAMH) algorithm that improved computational efficiency in sample-based imaging of electrical conductivity (EIT). That work used a linear approximation to the forward map in the first step of the algorithm. In this paper, we develop an alternative approximation for use in DAMH, namely a multilevel approximation developed from the hierarchy of coarse-scale models obtained by variational coarsening. This approach builds on two important strengths of robust multigrid solvers. First, the cost of a fine-scale solution of the forward map scales linearly with the degrees of freedom, and hence, it is provides better efficiency for algorithms performing sample-based inference. Second, the homogenization implicit in robust variational multigrid methods gives better solutions at coarse scales than homogenization by averaging of coefficients. We report results from a stylized example in electrical impedance imaging where data is a noisy and incomplete measurement of the Dirichlet-to-Neumann map.

A new efficient multilevel upscaling procedure for single-phase saturated flow in porous media is presented. While traditional approaches to this problem have focused on the computation of an upscaled hydraulic conductivity, here the coarse-scale model is created explicitly from the fine-scale model through the application of operator-induced variational coarsening. This technique, which originated with robust multigrid solvers, has been shown to accurately capture the influence of fine-scale heterogeneous structure over the complete hierarchy of coarse-scale models that it generates. Moreover, implicit in this hierarchy is the construction of interpolation operators that provide a natural and complete multiscale basis for the fine-scale problem. Thus, this new multilevel upscaling methodology is similar to the Multiscale Finite Element Method (MSFEM) and, indeed, it attains similar accuracy in computations of the fine-scale hydraulic head and coarse-scale normal flux on a variety of problems; yet it is an order of magnitude faster on the examples considered here.

We present a memory efficient parallel algorithm for the solution of tridiagonal linear systems of equations that are diagonally dominant on a very large number of processors. Our algorithm can be viewed as a parallel partitioning algorithm. We illustrate its performance using some examples. Based on this partitioning algorithm, we introduce a recursive version that has logarithmic communication complexity.

The main goal of this paper is to establish the convergence of mimetic discretizations of the first-order system that describes linear diffusion. Specifically, mimetic discretizations based on the support-operators methodology (SO) have been applied successfully in a number of application areas, including diffusion and electromagnetics. These discretizations have demonstrated excellent robustness, however, a rigorous convergence proof has been lacking. In this research, we prove convergence of the SO discretization for linear diffusion by first developing a connection of this mimetic discretization with Mixed Finite Element (MFE) methods. This connection facilitates the application of existing tools and error estimates from the finite element literature to establish convergence for the SO discretization. The convergence properties of the SO discretization are verified with numerical examples.

Minipermeameters are rapidly becoming a popular tool for collecting localized measurements of permeability in both laboratory and field studies. While one of the main advantages of minipermeameters is their ability to collect data on various support volumes, there have been only limited attempts to analyze their size and geometry. We define the support volume of minipermeameter measurements as a region containing 90% of the total gas flow, i.e., a region bounded by the 10% streamline. Using our new semi-analytical solutions for the Stokes' stream function we demonstrate that the support volume has the shape of a semi-toroid adjacent to the sample surface. Hence there is a blind spot directly below the minipermeameter, which is not probed by the measurement. We demonstrate that the support volume of the minipermeameter measurements decreases with the tip-seal's ratio (a ratio of the inner tip-seal radius to the outer tip-seal radius), while the size of the corresponding blind spot increases.

A broad class of discretizations of the diffusion operator is based on its first order form, allowing the rigorous enforcement of many desirable physical properties of the continuous model. In this research we investigate the development of multilevel solvers for the local or hybrid forms of these discretizations on logically rectangular quadrilateral meshes. In this case, the local elimination of flux leads to a system that contains both cell- and edge-based scalar unknowns. Based on this natural partitioning of the system we develop approximate reduced systems that reside on a single logically rectangular grid. Each such approximate reduced system, formed as an approximate Schur complement or as a variational product, are used as the first coarse-grid in a multigrid hierarchy or as a preconditioner for Krylov based methods.

The multiscale structure of heterogeneous porous media prevents a straightforward numerical treatment of the underlying mathematical flow models. In particular, fully resolved flow simulations are intractible and yet the fine-scale structure of a porous medium may significantly influence the coarse-scale properties of the solution (e.g., average flow rates). Consequently, homogenization or upscaling procedures are required to define approximate coarse-scale models suitable for efficient computation. Unfortunately, inherent in such a procedure is a compromise between its computational cost and the accuracy of the resulting coarse-scale solution. In general, most popular upscaling methods do not balance these competing demands. In this paper we highlight a new efficient, numerical method, which combines our recent work on multigrid homogenization (MGH) with the work of Dvořák (1994) to compute bounded estimates of the homogenized permeability for single phase saturated flows. Our approach is motivated by the observation that the coarse-scale influence of multiscale structures are captured automatically by robust variationally defined multigrid methods. The effectiveness of this new algorithm is demonstrated with numerical examples.

Avalon is a 140 processor Alpha/Linux Beowulf cluster constructed entirely from commodity personal computer technology and freely available software. Computational Physics simulations performed on Avalon resulted in the award of a 1998 Gordon Bell price/performance prize for significant achievement in parallel processing. Avalon ranked as the 113th fastest computer in the world on the November 1998 TOP500 list, obtaining a result of 47.8 Gigaflops on the parallel Linpack benchmark.

Avalon currently provides over 15,000 node-hours of production computing time per week, split among about 10 production users. Obtaining an equivalent amount of computing through Los Alamos institutional sources would cost a minimum of $45,000 per week. The machine also supports code development for another 50 users. The largest single simulation was performed in April and May 1998 using the SPaSM molecular dynamics code, which computed a total of 1.12 ×10¹⁶ floating point operations. This simulation is among the few scientific simulations to have ever involved more than 10 Petaflops of computation.

The price of hardware and final assembly labor for Avalon totalled $313,000 dollars. The monetary cost of the development and OS software used for the applications mentioned above was $0. Perhaps most extraordinary, all of the hardware and software maintenance on the machine is performed in the spare time of four people, averaging less than 10 man-hours of labor per week overall.

In mathematical models of flow through porous media, the coefficients typically exhibit severe variations in two or more significantly different length scales. Consequently, the numerical treatment of these problems relies on a homogenization or upscaling procedure to define an approximate coarse-scale problem that adequately captures the influence of the fine-scale structure. Inherent in such a procedure is a compromise between its computational cost and the accuracy of the resulting coarse-scale solution. Although techniques that balance the conflicting demands of accuracy and efficiency exist for a few specific classes of fine-scale structure (e.g., fine-scale periodic), this is not the case in general.

In this paper we propose a new, efficient, numerical approach for the homogenization of the permeability in models of single-phase saturated flow. Our approach is motivated by the observation that multiple length scales are captured automatically by robust multilevel iterative solvers, such as Dendy's black box multigrid. In particular, the operator-induced variational coarsening in black box multigrid produces coarse-grid operators that capture the essential coarse-scale influence of the medium's fine-scale structure. We derive an explicit local, cell-based, approximate expression for the symmetric, 2 × 2 homogenized permeability tensor that is defined implicitly by the black box coarse-grid operator. The effectiveness of this black box multigrid numerical homogenization method is demonstrated through numerical examples.

Certain classes of nodal methods and mixed-hybrid finite element methods lead to equivalent, robust, and accurate discretizations of second-order elliptic PDEs. However, widespread popularity of these discretizations has been hindered by the awkward linear systems which result. The present work overcomes this awkwardness and develops preconditioners which yield solution algorithms for these discretizations with an efficiency comparable to that of the multigrid method for standard discretizations. Our approach exploits the natural partitioning of the linear system obtained by the mixed-hybrid finite element method. By eliminating different subsets of unknowns, two Schur complements are obtained with known structure. Replacing key matrices in this structure by lumped approximations, we define three optimal preconditioners. Central to the optimal performance of these preconditioners is their sparsity structure which is compatible with standard finite difference discretizations and hence treated adequately with only a single multigrid cycle. In this paper we restrict the discussion to the two-dimensional case; these techniques are readily extended to three dimensions.

Nodal Methods have long been one of the most popular discretization techniques employed within the reactor physics community, while remaining conspicuously absent from the mainstream numerical analysis literature. A fundamental reason for this anomaly is that the physical arguments which were used to develop and enhance these methods seemed at odds with more rigorous discretization techniques. To facilitate communication between these distinct communities, a detailed chronological study of the lowest-order nodal methods for linear second order elliptic problems is presented. The presentation highlights the underlying motivation of these methods and formalizes many of their renowned properties. In addition, various equivalence relations within this family of discretizations are demonstrated, and equivalences with specific non-conforming and mixed-hybrid finite element methods (FEMs) are established. Rigorous error bounds and stability properties follow immediately from these latter equivalence relations, corroborating the results of a more rudimentary truncation error analysis.

An inherent difficulty of reactor simulation is that the coefficients in the mathematical model exhibit severe variations on two significantly different length scales. As in many other applications this is treated by defining an appropriate homogenization procedure which yields a simplified model with piecewise constant coefficients on a coarse scale suitable for efficient computation. Significant enhancements in accuracy are possible if the processes of homogenization and discretization are unified. We review the popular techniques that are based on this premise and rely on certain properties of the nodal discretization. In addition, we address the factors that contribute to their success in reactor modelling and deter their generalization outside of the reactor physics community. As an alternative to these highly specialized methods, we introduce a new multi-level homogenization technique which is readily applicable in a general setting, and is shown to have many important attributes.

Widespread acceptance of nodal methods has also been hindered by their use of nonstandard unknowns, as this results in stencils that appear awkward and incompatible with sophisticated iterative solution techniques. Specifically, equivalence with certain mixed-hybrid FEMs reveals that the nodal discretizations result in an indefinite system, which in two dimensions contains both cell-based and edge-based unknowns. Yet, inherent in this structure is a natural partitioning of the system which may be exploited to define a hierarchy of reduced systems (i.e. Schur complements) that are symmetric positive definite. Unfortunately, the reduced systems that involve unknowns of only one type, and hence seem most compatible with sophisticated iterative methods, suffer a loss of sparsity. However, the structure inherent in this hierarchy may be utilized in the construction of sparse approximate Schur complements for these systems. It is demonstrated that any one of these approximate operators, which are of either the standard 5 or 9-point family, may be utilized as an excellent preconditioner for conjugate gradient iterations. The efficiency of this approach is fully realized when the preconditioner is approximately inverted using only a single V- or W-cycle of a robust Black Box multigrid solver.

Measurements of the phase and modulation of amplitude-modulated light diffusely reflected by turbid media can be used to deduce absorption and scattering coefficients.

A time dependent diffusion model of photon migration in tissue is used to develop analytic expressions for the diffusely reflected pulse detected some distance from a delta function input. Particular attention is paid to the nature of the boundary between the tissue and the surrounding non-scattering medium, and it is shown that the pulse shape is relatively insensitive to the nature of this boundary. Monte Carlo simulation and experimental results are presented which confirm the accuracy of the diffusion model.

When a picosecond light pulse is incident on an optically turbid medium such as tissue, the temporal distribution of diffusely reflected and transmitted photons depends on the optical absorption and scattering properties of the medium. From diffusion theory it is possible to derive analytic expression for the pulse shape in terms of optical interaction coefficients of a homogeneous semi-infinite medium. Experimental tests of this model in tissue-simulating liquid phantoms of different geometries are presented here.

The increasing use of visible and near infrared light in therapeutic and diagnostic techniques has created a need to model its propagation in tissue. One of the fundamental objectives of such a model is the noninvasive evaluation of the optical properties of tissue. The focus of this thesis was the development of the diffusion approximation in the semi-infinite, slab, cylindrical and spherical geometries. This development required the derivation of approximate boundary conditions which included the zero, extrapolated and partial current boundary conditions. Calculations of the fluence and its related quantities arising from the extrapolated boundary condition were found to be in excellent agreement with the results of the more rigorous partial current boundary condition.

A preliminary evaluation of the validity of diffusion theory was performed by comparing its predictions to exact analytical calculations of the fluence in an infinite medium as well as Monte Carlo simulations of the reflectance and transmittance in 1-dimensional planar geometries. In all cases the agreement at late times was excellent. A practical test of the diffusion model was accomplished with the analysis of the reflectance data from a phantom of known optical properties in both the semi-infinite and slab geometries. The model perfomed well at low concentrations of added absorber, but a considerable discrepancy was found at the highest concentration. A systematic examination of the accuracy of the diffusion model as a function of the fundamental parameters is required to resolve this inconsistency.

Approximate expressions describing the equivalent information in the frequency domain were also developed for a semi-infinite medium. These expressions were then used to analyze the phase and modulation obtained from phantoms of known optical properties. Once again reasonable results were obtained at low concentrations of added absorber while a significant discrepancy arose at the highest concentration. The resolution of these descrepancies requires further investigation.

Since biological response to photodynamic therapy (PDT) depends on the light fluence distribution and photosensitizer concentration in the tissue, these two variables should ideally be measured noninvasively in individual cases. This can be reduced to determining the optical absorption and transport scattering coefficients of the tissue because, if these two parameters are known, the fluence distribution can be calculated and the photosensitizer concentration can be deduced from tis characteristic contribution to the absorption spectrum. The temporal spreading of a picosecond laser pulse as it propagates through tissue carries information about both the scattering and absorption properties of the tissue. A mathematical model is presented which allows derivation of the interaction coefficients from the pulse shape, and preliminary experiments are reported which demonstrate the potential of these techniques in PDT dosimetry. Equivalent information can be obtained in the frequency domain by using modulated light sources and measuring the phase and modulation of the detected light. Analytical expressions are developed for these observable quantities in terms of the optical interaction coefficients. Particular attention is drawn to the potential of low (less than 200MHz) frequency measurements as these can be made with relatively simple instrumentation.

The development of diagnostic and therapeutic photomedicine has generated a need to determine the optical properties of tissues in the U.V., visible and intrared regions of the spectrum. In this paper we will review the experimental techniques and resulting data on the optical prperties of mamalian tissues. These will include recent results from this laboratory as well as a summary of the work of other groups. Measurements are most abundant at around 630nm, the wavelength of greatest current interest for clinical photodynamic therapy. We will use these data as the reference values for examaning the wavelength-dependence of the optical properties.