TR92-43           R.A. Tapia and Yin Zhang

On the Quadratic Convergence of the Singular Newton's Method

The purpose of this essay is to describe a situation that we have found particularly exciting in our recent work in interior-point methods for linear programming. To our surprise, we have seen considerable theory developed concerning the superlinear convergence of singular Newton's methods. Hopefully, our comments will motivate further research in the general area of fast convergence for the singular Newton's method.

December 1992


TR92-42           C. Gonzaga and R.A. Tapia

On the Quadratic Convergence of the Simplified Mizuno-Todd-Ye Algorithm for Linear Programming

It is known that the Mizuno-Todd-Ye predictor-corrector primal-dual Newton interior-point method generates a duality gap sequence which converges quadratically to zero, and this is accomplished with an iteration complexity of O (sqrt(n) L). Very recently the present authors demonstrated that the iteration sequence generated by this method converges, and this convergence is to the analytical center of the solution set. In the current work we show that within a finite number of iterations the Newton corrector step can be replaced with a simplified Newton corrector step and the resulting algorithm maintains O (sqrt(n) L) iteration complexity, quadratic convergence of the duality gap sequence to zero, and convergence of the iteration sequence, however not necessarily to the analytical center. The simplified predictor-corrector algorithm requires only one linear solve per iteration in contrast to two linear solves per iteration required by the original predictor-corrector algorithm.

December 1992 (revised September 1994)

[SIAM J. Optim. 7 (1997), no. 1, 66-85.]


TR92-41           C. Gonzaga and R.A. Tapia

On the Convergence of the Mizuno-Todd-Ye Algorithm to the Analytic Center of the Solution Set

(Abstract not available)

December 1992 (revised September 1994)

[SIAM J. Optim. 7 (1997), no. 1, 47-65.]


TR92-40           A.S. El-Bakry and R.A. Tapia

On the Formulation of the Primal-Dual Newton Interior-Point Method for Nonlinear Programming

In this work we first study in detail the formulation of the primal-dual interior-point method for linear programming. We show that, contrary to popular belief, it cannot be viewed as the damped Newton's method applied to the Karush-Kuhn-Tucker conditioned for the logarithmic barrier function problem. Next we extend the formulation to general nonlinear programming, and then validate this extension by demonstrating that this algorithm can be implemented so that it is locally and Q-quadratically convergent under only the standard Newton's method assumptions. We discuss the globalization of the algorithm and include considerable numerical experimentation.

December 1992 (revised April 1995)


TR92-39           E.J. Cramer, P.D. Frank, G.R. Shubin, J.E. Dennis and R.M. Lewis

On Alternative Problem Formulations for Multidisciplinary Design Optimization

In this paper we introduce a perspective on multidisciplinary design optimization (MDO) problem formulation that provides a basis for choosing among existing formulations and suggests provocative, new ones. MDO problems offer a richer spectrum of possibilities for problem formulation than do single discipline design optimization problems, or multidisciplinary analysis problems. This is because the variables and the equations that characterize the MDO problem can be "partitioned" in some interesting ways between what we traditionally think of as the "analysis code(s)" and the "optimization code." An MDO approach can be characterized by what part of the overall computation is done in each code, how that computation is done, and what information is communicated between the codes.

December 1992


TR92-38           Philip T. Keenan

A Remark on Collocation and Upwinding in First Order Hyperbolic Systems

Keenan[2] defines and analyzes a new numerical method for coupled systems of nonlinear first order hyperbolic partial differential equations with one degenerate eigenvalue. That work extends in a certain direction the collocation method described by Luskin[3], which applies to systems with all the eigenvalues uniformly bounded away from zero. Luskin's method and Keenan's method both have direct application to the study of one dimensional fluid flow through pipelines. The pressure and velocity of an isothermal fluid in a pipeline can be described by a coupled pair of nonlinear first order hyperbolic partial differential equations. When thermal effects are important a third equation for temperature is added. While Luskin's method works well for the isothermal situation he discussed, it does not apply in certain common cases when thermal effects are modeled. The analysis of the new method shows how the difficulties that come from the application of standard collocation can be overcome by using upwinded piecewise constant functions for the degenerate component of the solution. Experiments indicate that this method is a substantial improvement over standard collocation.

A number of technical details obscure the analysis presented in Keenan [2], because that work treats the general nonlinear case. The present paper describes and analyzes the method in the context of a linear, constant coefficient system of equations. In this special case the proof simplifies considerably.

November 1992


TR92-37           Todd Arbogast and Mary F. Wheeler

A Characteristics-Mixed Finite Element Method for Advection Dominated Transport Problems

We define a new finite element method, called the characteristics-mixed method, for approximating the solution to an advection dominated transport problem. The method is based on a space-time variational form of the advection-diffusion equation. Our test functions are piecewise constant in space, and in time they approximately follow the characteristics of the advective (i.e., hyperbolic) part of the equation. Thus the scheme uses a characteristic approximation to handle advection in time. This is combined with a low order mixed finite element spatial approximation of the equation. Boundary conditions are incorporated in a natural and mass conservative fashion. The scheme is completely locally conservative; in fact, on the discrete level, fluid is transported along the approximate characteristics. A post-processing step is included in the scheme in which the approximation to the scaler unknown is improved by utilizing the approximate vector flux. This has the effect of improving the rate of convergence of the method. We show that it is optimally convergent to order one in time and at least suboptimally convergent to order 3/2 in space.

November 1992

[SIAM J. Numer. Anal. 32 (1995), no. 2, 404-424.]


TR92-36           R.E. Bixby, E. Andrew Boyd, Shireen S. Dadmehr and Ronni R. Indovina

The MIPLIB Mixed Integer Programming Library

(Abstract not available)

November 1992


TR92-35           W.W. Symes

A Consortium Proposal: "The Rice Inversion Project"

This document details a proposal for an industrially sponsored consortium for research in seismic inversion at Rice University. This consortium project will be directed by Professor William W. Symes in the Department of Computational and Applied Mathematics (formerly called The Department of Mathematical Sciences), George Brown School of Engineering. This project will develop novel approaches pioneered by Professor Symes to velocity and reflectivity estimation from waveform data, and will offer its sponsors both pilot software for state-of-the-art vector and parallel computing platforms, and a database of experience in waveform inversion. Participants will include graduate student assistants, postdoctoral research associates, and (whenever possible) short-and long-term visitors from the sponsoring organizations.

November 1992


TR92-34           Gang Bao and William W. Symes

Computation of Pseudo-Differential Operators

A simple algorithm is described for computing general pseudo-differential operator actions. Our approach is based on the asymptotic expansion of the symbol together with the Fast Fourier Transform (FFT). The idea is motivated by the characterization of pseudo-differential operator algebra. We show that the algorithm is efficient through analyzing its complexity. Some numerical experiments are also presented.

November 1992

[SIAM J. Sci. Comput. 17 (1996), no. 2, 416-429.]


TR92-33           M.L. Rosemblun

Automatic Differentiation: Overview and Application to Systems of Parameterized Nonlinear Equations

Automatic Differentiation is a computational technique that allows the evaluation of derivatives of functions defined by computer programs. Derivatives are calculated by applying the chain rule of differential calculus to the sequence of elementary computations involved in the program. In this work, an overview of the theory and implementation of automatic differentiation is presented, as well as a description of the available software.

An application of automatic differentiation in the context of solving systems of parameterized nonlinear equations is discussed. In this application, the "differentiated" functions are implementations of Newton's method and Broyden's method. The iterates generated by the algorithms are differentiated with respect to the parameters. The results show that whenever the sequence of iterates converges to a solution of the system, the corresponding sequence of derivatives (computed by automatic differentiation) also converges to the correct value. Additionally, we show that the "differentiated" algorithms can be successfully employed in the solution of parameter identification problems via the Black-Box method.

October 1992


TR92-32           C. Chiang, G. Raven and C. Dawson

How to Relate Monitoring Well and Aquifer Solute Concentrations

Field data show disparities between organic solute concentration in the aquifer and that in monitoring wells; an order of magnitude disparity has been recorded in some cases. Therefore, it is important to be able to relate concentrations between the two media for design of remediation systems. More significantly, to assess the impact of the leachate from a landfill on a downgradient drinking water well, it is important to correlate aquifer concentration with that in a drinking water well such that the true risk is not overly estimated.

A three-dimensional finite difference flow and transport model applied to demonstrate the disparity between aquifer and well concentrations is indeed well founded and can be quantified. The modeling shows that the concentration in the well is a function of the initial vertical concentration profile in the aquifer, the amount of flux from below the partially penetrated well, the degree of penetration, the soil lithology, and the amount of purged water before sampling. Based on these parameters, an approximate analytical solution is developed that agrees well with numerical solutions.

October 1992


TR92-31           A. Chilakapati and M.F. Wheeler

Convergence and Material Balance with MMOC-Galerkin

A Modified Method of Characteristics (MMOC) combined with Galerkin finite elements has often been used in the past to solve the advection dominated advection-diffusion equation that arises in miscible displacement, transport of soluble contaminants in groundwater and the bioremediation of contaminated aquifers. In this method the hyperbolic part of the equation is treated with characteristics and the remaining elliptic equation is solved with Galerkin finite elements. The right-hand-side in the later procedure is obtained through numerical integration. We demonstrate here that the error associated with this numerical integration is substantial enough to cause large material balance errors and poor convergence, even in the absence of overshoot and undershoot, typical of the Galerkin procedure. One can reduce the error by taking many quadrature points but at the cost of large CPU time. Finer discretization in space a and in time does not guarantee a better solution when the right-hand-side is not computed exactly . An exact integration scheme is implemented in 1-D. This is conservative and obtains theoretical convergence in the absence of overshoot and undershoot.

September 1992


TR92-30           M. El-Alem

A Robust Trust-Region Algorithm with a Non-Monotonic Penalty Parameter Scheme for Constrained Optimization

An algorithm for solving the problem of minimizing a non-linear function subject to equality constraints is introduced. This algorithm is a trust-region algorithm. In computing the trial step, a projected-Hessian technique is used that converts the trust-region subproblem to one similar to that of the unconstrained case. To force global convergence, the augmented Lagrangian is employed as a merit function.

One of the main advantages of this algorithm is the way that the penalty parameter is updated. We introduce an updating scheme that allows (for the first time, to the best of our knowledge) the penalty parameter to be decreased whenever it is warranted. The behavior of this penalty parameter is studied.

A convergence theory for this algorithm is presented. It is shown that this algorithm is globally convergent and that the globalization strategy will not disrupt fast local convergence. The local rate of convergence is also discussed. This theory is sufficiently general that it holds for any algorithm that generates steps whose normal components give at least a fraction of Cauchy decrease in the quadratic model of the constraints and uses Fletcher's exact penalty function as a merit function.

September 1992 (revised September 1993)

[SIAM J. Optim. 5 (1995), no. 2, 348-378.]


TR92-29           Catherine M. Samuelsen

The Dikin-Karmarkar Principle for Steepest Descent

Steepest feasible descent methods for inequality constrained optimization problems have commonly been plagued by short steps. The consequence of taking short steps is slow convergence to non-stationary points (zigzagging). In linear programming, both the projective algorithm of Karmarkar (1984) and its affined-variant, originally proposed by Dikin (1967), can be viewed as steepest feasible descent methods. However, both of these algorithms have been demonstrated to be effective and seem to have overcome the problem of short steps. These algorithms share a common norm. It is this choice of norm, in the context of steepest feasible descent, that we refer to as the Dikin-Karmarkar Principle.

This research develops mathematical theory to quantify the short step behavior of Euclidean norm steepest feasible descent methods and the avoidance of short steps for steepest feasible descent with respect to the Dikin-Karmarkar norm. While the theory is developed for linear programming problems with only nonnegativity constraints on the variables. Our numerical experimentation demonstrates that this behavior occurs for the more general linear program with equality constraints added. Our numerical results also suggest that taking longer steps is not sufficient to ensure the efficiency of a steepest feasible descent algorithm. The uniform way in which the Dikin-Karmarkar norm treats every boundary is important in obtaining a satisfactory convergence.

September 1992


TR92-28           J. E. Dennis, Jr., Mahmoud El-Alem, and Maria C. Maciel

A Global Convergence Theory for General Trust-Region-Based Algorithms for Equality Constrained Optimization

This work presents a global convergence theory for a broad class of trust-region algorithms for the smooth nonlinear programming problem with equality constraints. The main result generalizes Powell's 1975 result for unconstrained trust-region algorithms.

The trial step is characterized by very mild conditions on its normal and tangential components. The normal component must satisfy a fraction of Cauchy decrease condition on the quadratic model of the linearized constraints. The tangential component then must satisfy a fraction of Cauchy decrease condition on a quadratic model of the Lagrangian function in the translated tangent space of the constraints determined by the normal component. The Lagrange multipliers and the Hessians are assumed only to be bounded.

The other main characteristic of this class of algorithms is that the step is evaluated by using the augmented Lagrangian as a merit function and the penalty parameter is updated using the El-Alem scheme. The properties of the step together with the way that the penalty parameter is chosen are sufficient to establish global convergence.

As an example, an algorithm is presented which can be viewed as a generalization of the Steihaug-Toint dogleg algorithm for the unconstrained case. It is based on a quadratic programming algorithm that uses as a feasible point a step in the normal direction to the tangent space of the constraints and then does feasible conjugate reduced-gradient steps to solve the quadratic program. This algorithm should cope quite well with large problems for which effective preconditions are known.

September 1992 (revised August 1995)

[SIAM J. Optim. 7 (1997), no. 1, 177-207.]


TR92-27           C.N. Dawson, C.J. van Duijn and M.F. Wheeler

Characteristic-Galerkin Methods for Contaminant Transport with Non-Equilibrium Adsorption Kinetics

A procedure based on combining the method of characteristics with a Galerkin finite element method is analyzed for approximating reactive transport in groundwater. In particular, we consider equations modeling contaminant transport with nonlinear, non-equilibrium adsorption reactions. This phenomena gives rise to non-Lipschitz but monotone nonlinearities which complicate the analysis. A physical and mathematical description of the problem under consideration is given, then the numerical method is described and a priori error estimates are derived.

August 1992 (revised July 1993)

[SIAM J. Numer. Anal. 31 (1994), no. 4, 982-999.]


TR92-26           H. Song and G. Zhang

Potential Inversion of the Two Dimensional Plasma Wave Equation

A layer-stripping type method is developed for solving an inverse problem for the two dimensional plasma wave equation where the object is to find the potential given Cauchy data on a time-like surface. The key point of the method is that we use one way wave approximation equations instead of the full wave equation for extrapolation of the wave field to avoid instability in numerical computation. Numerical experiments are performed to examine the effectiveness of the method. Numerical results show that this method works well for synthetic data.

August 1992

[Inverse Problems 10 (1994), no. 2, 467-484.]


TR92-25           H. Song and G. Zhang

Wavefield Splitting and Extrapolation of the Two Dimensional Plasma Wave Equation

Wave splitting is an important technique in a great variety of studies concerning wave propagation. It is a class of methods of splitting the wave equation under consideration into one way wave equations, or physically splitting the wave field into two components with opposite propagation directions. In this paper, a new approach for wave splitting is presented to get a coupled one way wave system for the two dimensional plasma wave equation by using the theory and techniques of pseudo-differential operators. The coupled system and the original equation are equivalent in the sense that they are the same for the singularities propagating in non-glancing directions. A localized approximation of the nonlocal system is given and energy estimates of some related wavefield extrapolation problems, corresponding to the migration problem and wavefield downward continuation problem in exploration geophysics respectively, are obtained. An application of the system to an inverse potential problem is outlined.

August 1992

[In Mathematical and Numerical Aspects of Wave Propagation (Mandelieu-La Napoule, 1995), 328-337, SIAM, Philadelphia, PA, 1995.]


TR92-24           C.N. Dawson and T.F. Dupont

Several Procedures for Operator-Based Averaging for Elliptic Equations

Numerical procedures are discussed for constructing averaged coefficients for elliptic differential operators. These procedures are intended for problems where the coefficients vary on a scale finer than can be resolved by a reasonable computational grid. Numerical methods for calculating locally averaged coefficients using mixed and Galerkin finite elements are presented. These methods involve solving local elliptic problems either to determine a pseudo-coefficient or as part of the overall solution procedure. The local problems are independent and can be solved in parallel. The procedures are formulated and numerical results demonstrating their performance are presented.

August 1992


TR92-23           W.W. Symes

DSO Velocity Inversion: A "Gas Cloud" Synthetic Example

A low velocity near surface anomaly delays the seismic wavefield and focuses wavefronts passing through it. Such anomalies occur as a result of localized gas seepages ("gas clouds") in the sedimentary column, for example. Use of a migration velocity model not containing the anomaly will produce images containing false synclinal structures. Application of DSO velocity inversion to a simple example of this type images the velocity anomaly and removes the false structure.

August 1992


TR92-22           W.W. Symes

The Plane-Wave Detection Problem

The plane-wave detection problem is: to estimate the incidence angle and waveform of a transient plane traveling wave, from samples recorded at a linear array of receivers. This simple problem shares several important mathematical features with other inverse problems of wave propagation, and is of interest in its own right as a model problem in ocean acoustic signal analysis. Straightforward formulation as a nonlinear least squares problem yields a nonconvex objective for which the minima are not stably dependent on the data. In contrast, an infeasible point formulation, in which the signal at each receiver is explained to some extent independently, proves to yield a smooth convex optimization problem with stable optima. Numerical experiments illustrate the theoretical results about the infeasible point approach, differential semblance optimization.

August 1992 (revised April 1993)

[Inverse Problems 10 (1994), no. 6, 1361-1391.]


TR92-21           J.E. Dennis and R.M. Lewis

Nonlinear Programming and Domain Decomposition for Partial Differential Equations

(Abstract not available)

August 1992


TR92-20           H. Yabe

Variations of Structured Broyden Families for Nonlinear Least Squares Problems

In this paper, we consider structured quasi-Newton methods for finding a local solution to nonlinear least squares problems. This paper is concerned with the line search globalization method. Recently, factorized versions of the structured quasi-Newton methods have been studied by Sheng and Zou, and Yabe and Takahashi in order to obtain a descent search direction for the objective function. In this paper we first generalize the update of Sheng and Zou, and propose a new factorized family corresponding to the Broyden family (SZ-Broyden family). Second, we suggest a relationship between the structured quasi-Newton updates and the factorized versions. We use this relationship to show that the factorized Broyden family proposed by Yabe and Yamaki corresponds to the Engels and Martinez family, and we further obtain a new structured quasi-Newton update which corresponds to the SZ-Broyden family in the sense of this relation. Finally, we apply sizing techniques to these methods and present some numerical experiments.

June 1992


TR92-19           Y. Zhang and A.S. El-Bakry

A Modified Predictor-Corrector Algorithm for Locating Weighted Centers in Linear Programming

In certain applications of linear programming, some particular solutions called the weighted centers of the solution set are often desired, giving rise to the need of algorithms capable of locating such centers. In this note, we modify the Mizuno-Todd-Ye predictor-corrector algorithm so that the modified algorithm is guaranteed to converge to the weighted center for given weights. The key idea is to ensure that iterates remain in a sequence of shrinking neighborhoods of the weighted central path. The modified algorithm also possesses polynomiality and superlinear convergence.

June 1992

[J. Optim. Theory Appl. 80 (1994), no. 2, 319-331.]


TR92-18           W.W. Symes

A Differential Semblance Criterion for Inversion of Multioffset Seismic Reflection Data

Mean-square error leading to least-squares inversion of multioffset reflection seismograms is insensitive to velocity trend information except in the immediate vicinity of a kinematically correct model. In contrast, differential semblance retains sensitivity to velocity trend changes over a wide range of models. The differential semblance criterion combines mean-square error with the mean-square differences of inverted models from datasets at neighboring shot positions (or offsets, or slownesses, ...). Differential semblance compares model estimates at nearby acquisition parameters which are similar even when the model velocity trends are incorrect. Because the method inverts the data, so that the estimated model amplitudes are meaningful, simple differences between the (unstacked) model estimates give a reliable measure of velocity error. A mathematical investigation indicates that the differential semblance criterion is smooth and convex over a large range of velocity models. Numerical simulation using synthetic data sets verifies this contention.

May 1992


TR92-17           R. Grundy, C.J. van Duijn and C.N. Dawson

Asymptotic Profiles with Finite Mass in One-Dimensional Contaminant Transport Through Porous Media: The Fast Reaction

(Abstract not available)

April 1992

[Quart. J. Mech. Appl. Math. 47 (1994), no. 1, 69-106.]


TR92-16           William W. Symes and Michel Kern

Velocity Inversion by Differential Semblance Optimization for 2D Common Source Data

Differential Semblance Optimization ("DSO") is a variant of data-fitting (least-squares) inversion of reflection seismograms. The misfit functions used in other implementations of least-squares inversion (for example Tarantola 1986, Kolb et al. 1986) exhibit highly nonconvex dependence on velocity trends. Therefore least-squares inversion appears to require the use of relatively costly global optimization algorithms such as Monte-Carlo minimization, simulated annealing, or genetic algorithms (Sen and Stoffa 1991, Scales et al. 1991, Tarantola 1991). Also local sensitivity analysis (eg. singular value decomposition of the Hessian) does not describe the resolution limits inherent in such misfit functions - they are too nonlinear. In contrast, the DSO misfit function is smooth and convex over a wide range of velocity models, and so can be minimized satisfactorily via efficient local (Newton-like) algorithms. Moreover the quadratic model of the DSO misfit function at the optimum is descriptive of its local behaviour, and so may be used for studies of velocity resolution.

In previous work (Symes and Carazzone 1991 and references cited there) we have presented DSO inversion for plane-wave data sets and layered media, including application to field data. In this paper, we describe the DSO misfit function for 2D shot gather inversion of laterally heterogeneous models, and an algorithm for minimizing it, and present a simple, preliminary example of shot-gather velocity inversion.

April 1992


TR92-15           M.F. Wheeler, K.R. Roberson, and Ashokkumar Chilakapati

Three-Dimensional Bioremediation Modeling in Heterogeneous Porous Media

Recently Rice University and Pacific Northwest Laboratory (PNL) have begun a collaborative research effort that involves laboratory, field, and simulation work directed toward validating remediation strategies, including both natural and in situ bioremediation at U.S. Department of Energy (DOE) sites such as Hanford. Because of chemical, biological, geologic and physical complexities of modeling these DOE sites, one of the major simulation goals of the project is to formulate and implement accurate and efficient (parallel) algorithms for modeling multiphase/multicomponent flow and reactive transport.

In this paper we first describe the physical problem that needs to be modeled. Because of the emergence of concurrent supercomputing, we propose accurate numerical algorithms that are based on operator-splitting in time and domain decomposition iterative techniques. In particular, three characteristic finite element methods and two operator-splitting algorithms are described for advection-diffusion-reaction problems. Three-dimensional bioremediation modeling results in a heterogeneous saturated porous medium are presented.

April 1992


TR92-14           Brian D. Wood and Clint N. Dawson

Effects of Lag and Maximum Growth in Contaminant Transport and Biodegradation Modeling

The effects of time lag and maximum microbial growth on biodegradation in contaminant transport are discussed. A mathematical model is formulated that accounts for these effects, and a numerical case study is presented that demonstrates how lag influences biodegradation.

April 1992


TR92-13           Philip T. Keenan

Mixed Methods on Quadrilaterals and Hexahedra

We describe a new family of discrete spaces suitable for use with mixed methods on general quadrilateral and hexahedral elements. The new spaces are natural in the sense of differential geometry, so all the usual mixed method theory, including the hybrid formulation, carries over to these new elements with proofs unchanged. Because transforming general quadrilaterals into squares introduces nonlinearity and because mixed methods involve the divergence operator, the new spaces are more complicated than either the corresponding Raviart-Thomas spaces for rectangles or corresponding finite element spaces for quadrilaterals. These new elements may be useful in topologically regular grids, where initially rectangular grids are deformed to match features of the physical region.

April 1992 (revised May 1992)

[Numer. Algorithms 7 (1994), no. 2-4, 269-293.]


TR92-12           Clint N. Dawson and Mary F. Wheeler

Time-Splitting Methods for Advection-Diffusion-Reaction Equations Arising in Contaminant Transport

Two time-splitting methods for advection-diffusion-reaction problems are discussed and analyzed. Numerical results for the first approach applied to bioremediation of contaminants in groundwater are presented.

March 1992

[In ICIAM 91 (Washington, DC, 1991), 71-82, SIAM, Philadelphia, PA, 1992.]


TR92-11           T. Arbogast

Gravitational Forces in Dual-Porosity Systems II: Computational Validation of the Homogonized Model

Three models are considered for single component, single phase flow in naturally fractured porous media. The microscopic model holds on the Darcy scale, and it is considered to govern the system. The macroscopic, dual-porosity model was derived in Part I of their work from the microscopic model by two-scale mathematical homogenization. In this paper, we show that the dual-porosity model predicts well the behavior of the microscopic model by comparing their computed solutions in certain reasonable test cases. Homogenization gives a complex formula for a key parameter in the dual-porosity model; herein a simple approximation to this formula is presented. The third model considered is a single-porosity model with averaged parameters. It is shown that this type of model cannot predict the behavior of the microscopic flow.

March 1992


TR92-10           Todd Arbogast

Gravitational Forces in Dual-Porosity Systems I. Model Derivation by Homogenization

We consider the problem of modeling flow through naturally fractured porous media. In this type of media, various physical phenomena occur on disparate length scales, so it is difficult to properly average their effects. In particular, gravitational forces pose special problems. In this paper we develop a general understanding of how to incorporate gravitational forces into the dual-porosity concept. We accomplish this through the mathematical technique of formal two-scale homogenization. This technique enables us to average the single-porosity, Darcy equations that govern the flow on the finest (fracture thickness) scale. The resulting homogenized equations are of dual-porosity type. We consider three flow situations, the flow of a single component in a single phase, the flow of two fluid components in two completely immiscible phases, and the completely miscible flow of two components.

March 1992


TR92-09           Virginia Torczon

PDS: Direct Search Methods for Unconstrained Optimization on Either Sequential or Parallel Machines

PDS is a collection of Fortran subroutines for solving unconstrained nonlinear optimization problems using direct search methods. The software is written so that execution on sequential machines is straightforward while execution on Intel distributed memory machines, such as the iPSC/2, the iPSC/860 or the Touchstone Delta, can be accomplished simply by including a few well-defined routines containing calls to Intel-specific Fortran library calls. Those interested in using the algorithm on other distributed memory machines, even for something as simple as a network of workstations or personal computers, need only modify these few subroutines to handle the global communication requirements. Furthermore, since the parallelism is clearly defined at the "do-loop" level, it is a simple matter to insert compiler directives that allow for execution on shared memory parallel machines. Included here is an example of such directives, contained in comment statements, for execution on a Sequent Symmetry S81.

PDS encompasses an entire class of general-purpose optimization methods that require little of the user other than a (scalar) subroutine to evaluate the function (though the algorithm is flexible enough to accommodate subroutines that evaluate the function in parallel) and even less of the problem to be solved since direct search methods presumes only that the function is continuous. Thus, these methods are particularly effective on parameter estimation problems involving a relatively small number of parameters. They are also very interesting as parallel algorithms because they are perfectly scalable: they can use any number of processors regardless of the dimension of the problem to be solved and, in fact, tend to perform better as more processors are added.

March 1992


TR92-08           E. Andrew Boyd

Solving Integer Programs with Enumeration Cuts

(Abstract not available)

March 1992


TR92-07           Mary F. Wheeler and John R. Whiteman

Superconvergence of Recovered Gradients of Discrete Time/Piecewise Linear Galerkin Approximations for Linear and Nonlinear Parabolic Problems

Superconvergent error estimates in l2(H¹) and linfinity(H¹) norms are derived for recovered gradients of finite difference in time/piecewise linear Galerkin approximations in space for linear and quasi-nonlinear parabolic problems in two space dimensions. The analysis extends previous results for elliptic problems to the parabolic context, and covers problems in regions with non-smooth boundaries under certain assumptions on the regularity of the solutions.

March 1992

[Numer. Methods Partial Differential Equations 10 (1994), no. 3, 271-294.]


TR92-06           Todd Arbogast, Ashokkumar Chilakapati and Mary F. Wheeler

A Characteristic-Mixed Method for Contaminant Transport and Miscible Displacement

Recently, Arbogast and Wheeler have formulated and analyzed a modified method of characteristics-mixed method for approximating solutions to convection-diffusion equations. This scheme is theoretically mass conservative over each grid cell; it is approximately so in implementations. We consider application of this procedure to contaminant transport and to miscible displacement with unfavorable mobility ratio. Results in one, two, and three space dimensions are discussed.

February 1992

[In Computational Methods in Water Resources, IX, Vol.1 (Denver, CO, 1992), 77-84, Comput. Mech., Southampton, 1992.]


TR92-05           Todd Arbogast

Simplified Dual-Porosity Model for Two-Phase Flow

A model for two-phase, incompressible, immiscible fluid flow in a highly fractured porous medium is derived as a simplification of a much more detailed dual-porosity model. This simplified model has a nonlinear matrix-fracture interaction, and it is more general than similar existing "transfer function" models. It is computationally less complex than the detailed model, and simulation results are presented which assess any loss in accuracy. It is shown that the new model approximates capillary effects quite well, and better than similar existing models.

February 1992

[In Computational Methods in Water Resources, IX, Vol. 2 (Denver, CO, 1992), 419-426, Comput. Mech., Southampton, 1992.]


TR92-04           Todd Arbogast

The Existence of Weak Solutions to Single Porosity and Simple Dual-Porosity Models of Two-Phase Incompressible Flow

It is shown that there exists a weak solution to a degenerate and singular elliptic-parabolic partial integro-differential system of equations. These equations model two-phase incompressible flow of immiscible fluids in either an ordinary porous medium or in a naturally fractured porous medium. The full model is of dual-porosity type, though the single porosity case is covered by setting the matrix-to-fracture flow terms to zero. This matrix-to-fracture flow is modeled simply in terms of fracture quantities; that is, no distinct matrix equations arise. The equations are considered in global pressure formulation that is justified by appealing to a physical relation between the degeneracy of the wetting fluid's mobility and the singularity of the capillary pressure function. In this formulation, the elliptic and parabolic parts of the system are separated; hence, it is natural to consider various boundary conditions, including mixed nonhomogeneous, saturation dependent ones of the first three types. A weak solution is obtained as a limit of solutions to discrete time problems. The proof makes no use of the corresponding regularized system. The hypotheses required for various earlier results on single-porosity systems are weakened so that only physically relevant assumptions are made. In particular, the results cover the cases of a singular capillary pressure function, a pure Neumann boundary condition, and an arbitrary initial condition.

February 1992

[Nonlinear Anal. 19 (1992), no. 11, 1009-1031.]


TR92-03           W.W. Symes and M. Kern

Inversion of Reflection Seismograms by Differential Semblance Analysis: Algorithm Structure and Synthetic Examples

Seismograms predicted from acoustic or elastic earth models depend very nonlinearly on the long wavelength components of velocity. This sensitive dependence demands the use of special variational principles in waveform-based inversion algorithms. The differential semblance variational principle is well-suited to velocity inversion by gradient methods, since its objective function is smooth and convex over a large range of velocity models. An extension of the adjoint state technique yields an accurate estimate of the differential semblance gradient. Nonlinear conjugate gradient iteration is quite successful in locating the global differential semblance minimum, which is near the ordinary least squares global minimum when coherent data noise is small. Several examples based on the 2D primaries-only acoustic model illustrate features of the method and its performance.

July 1992


TR92-02           Guangye Li and Irvin J. Lustig

An Implementation of a Parallel Primal-Dual Interior Point Method for Multicommodity Flow Problems

An implementation of the primal-dual predictor-corrector interior point method is specialized to solve linear multicommodity flow problems. The block structure of the constraint matrix is exploited via parallel computation. The bundling constraints require the Cholesky factorization of a dense matrix, where a method that exploits parallelism for the dense Cholesky factorization is used. The resulting implementation is 65 to 90 percent efficient, depending on the problem instance. For a problem with K commodities, an approximate speedup for the interior point method of 0.8K is realized.

January 1992


TR92-01           Guangye Li

A Diagonal-Secant Update Technique for Sparse Unconstrained Optimization

This paper presents a diagonal-secant modification of the successive element correction method, a finite-difference based method, for sparse unconstrained optimization. This new method uses the gradient values more efficiently in forming the approximate Hessian than the successive element correction method. It is shown that the new method has at least the same local convergence rates as the successive element correction method for general problems and that it has better q-convergence and f-convergence rates than the successive element correction method for problems with band structures. The numerical results show that the new method may be competitive with most of the existing methods for some problems.

January 1992


Go to 1991 reports

Go to 1993 reports