University of Southern California Dissertations and Theses

Algorithms for stochastic Galerkin projections: solvers, basis adaptation and multiscale modeling and reduction

(USC Thesis Other)

Algorithms for stochastic Galerkin projections: solvers, basis adaptation and multiscale modeling and reduction

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Algorithms for stochastic Galerkin projections: Solvers, basis adaptation and multiscale modeling and reduction by Ramakrishna Tipireddy A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (CIVIL ENGINEERING) August 2013 Copyright 2013 Ramakrishna Tipireddy Dedication To my mother, my first teacher. ii Acknowledgements Foremost, I would like to express my sincere gratitude to my advisor Prof. Roger Ghanem for his continuous support, patience, insightful suggestions and critical comments, enthusiasm and continuous push for betterment of the work throughout the course of my PhD. I would like to thank him for introducing me to wide areas of research and providing many opportunities to col- laborate with the experts in various research fields and present my work in numerous international conferences. His high academic standards will always be a source of inspiration for me. I would like to thank Dr. Eric Phipps for mentoring me during my internship at Sandia Na- tional Labs, for serving on my qualifying committee and providing valuable support in technical work and technical writing during the course of this work. I would like to thank my collaborators, Prof. Somnath Ghosh and Dr. Daniel Paquet for providing experimental data and their support during this work. I would like to thank my dissertation committee members, Prof. Sami Masri, Prof. Aiichiro Nakano and Prof. Robert Lucas for their valuable comments and suggestions for betterment of this work. I would like to thank Prof. Erik Johnson for serving on my qualifying committee and providing valuable suggestions regarding my presentation skills. I gratefully acknowledge the financial support received from US Department of Energy, Na- tional Science foundation during the course of this work. I would like to thank Viterbi school of engineering for offering me Viterbi fellowship. iii I would like to extend my sincere thanks to my former and current colleagues, friends who supported directly or indirectly during this work. I would like thank Dr. Maarten Arnst, Dr. Son- joy Das and Dr. Johann Guilleminot for their support and guidance in the earlier phase of my PhD. My sincere thanks to my former colleagues Dr. Arash Noshadravan, Dr. Maud Comboul, Dr. Sara Abedi, Daniel Lakeland, Dr. Evangelia Kalligiannaki for their input during many dis- cussions I had. I would also like to thank my other past and current fellow group members, Dr. Hadi Meidani, Dr. Yongbo Peng, Dr. Bedrick Sousdik, Iman Yadegaran, Vahid Keshavarzzadeh, Hamed Haddadzadegan, Charanraj Thimmisetty, Prasanth Koganti and Matthew Sexton for all the discussions and sharing thoughts on interesting problems. I would like to thank my friends, Dr. Chandrika Prakash Vyasarayani, Dr. Nagesh Janapala, Dr. Karthikeyan Chockalingam, Sukhpreet Sandhu, Dr. S. Kumar and Dr. Viswanath Chinthapenta for their constant support and encouragement. I would like to thank Prof. C. S. Manohar, who taught me various subjects, guided my research work at the Indian Institute of Science. He has been a source of inspiration and he motivated and encouraged me to pursue higher studies. I am grateful to all my teachers who taught me and inspired me throughout my student life. Last but not the least, I would like to thank my family: my father Chalama Reddy, my brother Harikrishna Reddy for their unconditional love and sacrifices without which this work would not have been completed. I would like to thank my uncle Kesavula Reddy, who has been influential in my early stages of education. iv Table of Contents Dedication ii Acknowledgements iv List of Tables viii List of Figures ix Abstract xiii Chapter 1 Introduction 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Stochastic upscaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Challenges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Chapter 2 Stochastic Galerkin methods 8 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Stochastic partial differential equations . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Input random field model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.1 Karhunen-Lo` eve expansion . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.2 Polynomial chaos expansion . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 Stochastic Galerkin method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Chapter 3 Iterative solvers and preconditioners 16 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Solution methods for stochastic Galerkin systems . . . . . . . . . . . . . . . . . 16 3.2.1 Jacobi mean algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2.2 Gauss-Seidel mean iterative method . . . . . . . . . . . . . . . . . . . . 18 3.2.3 Krylov based iterative methods with matrix-free operations . . . . . . . . 21 3.3 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.4 Sparse grid collocation method . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.5 Numerical illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Chapter 4 Basis adaptation - dimension reduction 39 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 v 4.2 Polynomial Chaos Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2.1 Change of Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.2.2 Reduction via Projection . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3 Special Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.3.1 Case where Gaussian components ofq are known . . . . . . . . . . . . . 45 4.3.2 Case where quadratic components ofq are known . . . . . . . . . . . . 46 4.4 Application to Equations with Random Coefficients . . . . . . . . . . . . . . . . 48 4.4.1 Non-Intrusive Implementation . . . . . . . . . . . . . . . . . . . . . . . 48 4.4.2 Intrusive Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.4.3 Adapted Representation of Model Parameters . . . . . . . . . . . . . . . 52 4.5 Numerical Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.5.1 Algebraic Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.5.2 Linear Elasticity with Nonlinear Quantities of Interest . . . . . . . . . . 55 4.5.3 Random eigenvalue problem . . . . . . . . . . . . . . . . . . . . . . . . 59 Chapter 5 Stochastic upscaling 66 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.2.1 Navier-Stokes equations at fine-scale . . . . . . . . . . . . . . . . . . . 68 5.2.2 Darcy-Brinkman equations at coarse-scale . . . . . . . . . . . . . . . . . 69 5.3 Stochastic upscaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.3.1 Permeability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.3.2 Thermal conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.4 Basis adaptation to find QoI . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.4.1 Linear basis adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.4.2 Quadratic basis adaptation . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.4.3 Basis adaptation with KL expansion . . . . . . . . . . . . . . . . . . . . 74 5.5 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 Chapter 6 High resolution micrograph synthesis 81 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 6.3 Parametric texture model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.3.1 Steerable pyramid decomposition . . . . . . . . . . . . . . . . . . . . . 90 6.3.2 Parameter estimation sketch . . . . . . . . . . . . . . . . . . . . . . . . 92 6.3.3 Micrograph synthesis procedure . . . . . . . . . . . . . . . . . . . . . . 95 6.4 Density-based Monte Carlo filter . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.5 Numerical illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.6 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 Bibliography 109 vi List of Tables 3.1 Scaled solution time (# of iterations) vs stochastic dimension for random field with uniform random variables, diffusion equation, PC order = 4 and=0.1 . . . 32 3.2 Scaled solution time (# of iterations) vs order of polynomial chaos for random field with uniform random variables, diffusion equation, Stoch. dim=4, = 0:1 . 33 3.3 Scaled solution time (# of iterations) vs standard deviation () of input random field, for random field with uniform random variables, diffusion equation, Stoch dim = 4 and PC order = 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.4 Scaled solution time (# of iterations) vs stochastic dimension for random field with uniform random variables, advection-diffusion equation, PC order = 4, and =0.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.5 Scaled solution time (# of iterations) vs order of polynomial chaos for random field with uniform random variables, advection-diffusion equation, Stoch. dim=4, = 0:1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.6 Scaled solution time (# of iterations) vs standard deviation () for random field with uniform random variables, advection-diffusion equation, Stoch dim = 4, PC order = 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.7 Scaled solution time (# of iterations) vs stochastic dimension for log-normal ran- dom field, diffusion equation, PC order = 4, and=0.1 . . . . . . . . . . . . . . 36 3.8 Scaled solution time (# of iterations) vs order of polynomial chaos for log-normal random field, diffusion equation, Stoch. dim=4, = 0:1 . . . . . . . . . . . . . 36 3.9 Scaled solution time (# of iterations) vs standard deviation () for log-normal random field, diffusion equation, Stoch dim = 4 and PC order = 4 . . . . . . . . . 37 3.10 Sacled solution time (# of iterations) vs stochastic dimension for log-normal ran- dom field, advection-diffusion equation, PC order = 4 and=0.1 . . . . . . . . . 37 3.11 Scaled solution time (# of iterations) vs order of polynomial chaos for log-normal random field, advection-diffusion equation, Stoch. dim=4, = 0:1 . . . . . . . . 38 vii 3.12 Scaled solution time (# of iterations) vs standard deviation () for log-normal random field, advection-diffusion equation, Stoch dim = 4 and PC order = 4 . . . 38 viii List of Figures 1.1 Fine scale domain modeled with an array of solid discs and coarse scale domain modeled as porous media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.1 Scaled Galerkin solution time vs solution error for varying polynomial chaos order and spatial mesh for the diffusion problem with uniform random field . . 30 3.2 Relative residual norm vs iteration countto solve stochastic diffusion equation . 31 3.3 Relative residual norm vs iteration count to solve stochastic advection-diffusion equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.1 Probability density function of h(). Results using a full 10-dimensional so- lution and the reduced approximations using prior knowledge on the Gaussian content and quadratic content of QoI. . . . . . . . . . . . . . . . . . . . . . . . 54 4.2 Elastic domain with random Young’s modulus and subjected to uniform internal pressure. PointA (0.2977,0.1417) is highlighted. . . . . . . . . . . . . . . . . 56 4.3 Probability density function of x–displacement at point A. Results using a 1) full 10-dimensional solution and the approximations using the 2) dominant component in the Karhunen-Loeve expansion of the stochastic process associ- ated with the solution, 3) the 1-d approximation associated with Gaussian prior knowledge, and 4) sequence of 10 1-d approximations associated with quadratic prior knowledge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.4 Coefficients of the PC expansion of the x-displacement at point A. Results using a full 10-dimensional solution and the reduced approximation using prior knowledge on the quadratic content of QoI. . . . . . . . . . . . . . . . . . . . 58 4.5 Probability density function of the maximum vonMises stress over the spatial domain of the problem. Results using 1) full 10-dimensional solution and the approximations using the 2) 1-d approximation associated with Gaussian prior knowledge, and 3) sequence of 10 1-d approximations associated with quadratic prior knowledge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 ix 4.6 The first six coefficients of the elasticity modulus with respect to a polynomial chaos expansion in the original gaussian variables. . . . . . . . . . . . . . . . 60 4.7 The first six coefficients of the elasticity modulus with respect to a polynomial chaos expansion in the one-dimensional gaussian variable 1 . . . . . . . . . . . 60 4.8 2-d cantilever beam meshed with quad4 elements, is fixed on left end and eigen- value analysis is performed . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.9 Wind turbine blade meshed with quad4 elements, is fixed in all 6 dof on one end and eigenvalue analysis is performed . . . . . . . . . . . . . . . . . . . . . . . 63 4.10 pdf of eigenvalues of 2-d cantilever beam in 6d and in 1d . . . . . . . . . . 64 4.11 pdf of eigenvalues of wind turbine blade in 6d and in 1d . . . . . . . . . . . 65 4.12 Polynomial chaos coefficients of of (l) j in 6d and in 1d . . . . . . . . . . . 65 5.1 Fine scale domain modeled with an array of solid discs and coarse scale domain modeled as porous media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2 Fluid and solid regions meshed with 2-d quad4 elements . . . . . . . . . . . . 76 5.3 Realizations of the thermal conductivity in the circular discs . . . . . . . . . . 76 5.4 Mean and std. dev of velocity in x-direction in 10d and order 3 . . . . . . . . 77 5.5 Mean and std. dev of velocity in y-direction in 10d and order 3 . . . . . . . . 77 5.6 Mean and std. dev of pressure in 10d and order 3 . . . . . . . . . . . . . . . 77 5.7 Mean and std. dev of temperature in 10d and order 3 . . . . . . . . . . . . . 77 5.8 Mean and std. dev. of permeability componentk 11 using volume averaging . . 78 5.9 Eigenevalues obianed from KL expansion of thermeal conductivity at finescale and the eigen values obtained from KL analysis of Gaussian part of permeability at coarse scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.10 Permeability computed at points 3 and 20, using linear and quadratic basis adap- tation methods. Here, the basis adapted for each spatial point (QoI) separately . 79 5.11 Permeability computed at points 3 and 20 using basis adaptation method, where new random variables are obtained from Karhunen-Loeve analysis of the Gaussian part of the random filed. Here, basis is adapted for all spatial points (QoI) simultaneously. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 x 5.12 Thermal conductivity computed at points 20 using basis adaptation method, where new random variables are obtained from Karhunen-Loeve analysis of the Gaussian part of the random filed. Here, basis is adapted for all spatial points (QoI) simultaneously. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 6.1 Low resolution micrograph Low resolution, low magnification digital image of cast aluminum alloy W319, for which high resolution images are available at points A and B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 6.2 High resolution digital image of cast aluminum alloy W319 at A and B . . . 86 6.3 Low resolution digital image of cast aluminum alloy W319 at A, B and C . 87 6.4 Steerable pyramid diagram Block diagram for the steerable pyramid Portilla and Simoncelli (2000), Simoncelli and Freeman (1995). Input image is initially split into high-pass and low-pass bands. Low-pass band is further decomposed into a set of oriented subbands and low-pass residual band. The recursive de- composition is achieved by inserting the diagram in dotted lines into the solid circle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.5 Real parts of subbands Real parts of steerable pyramid orientation bands(4 ori- entations and 4 scales) and low-pass subband of image at location A. High-Pass subband is not shown here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.6 Block diagram of the recursive micrograph synthesis procedure . . . . . . 97 6.7 Block diagram of synthesis procedure with particle filter. Block diagram of the recursive micrograph synthesis procedure along with the density-based Monte Carlo filter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.8 Micrograph at A synthesized without particle filter. Micrographs from left to right show respectively, the experimental high resolution micrograph, synthe- sized micrographs without particle filter which are statistically equivalent with high resolution micrograph at location A. . . . . . . . . . . . . . . . . . . . . 102 6.9 Micrograph at A synthesized with particle filter. From 6.9(a) to 6.9(c), experi- mental high resolution micrograph, experimental low resolution micrograph and high resolution micrograph synthesized with texture model and particle filter at location A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.10 Micrograph at B synthesized without particle filter. Micrographs from left to right show respectively, the experimental high resolution micrograph, synthe- sized micrographs without particle filter which are statistically equivalent with high resolution micrograph at location B. . . . . . . . . . . . . . . . . . . . . . 103 xi 6.11 Micrograph at B synthesized with particle filter. From6.11(a) to 6.11(c), experi- mental high resolution micrograph, experimental low resolution micrograph and high resolution micrograph synthesized with texture model and particle filter at location B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6.12 Micrograph at C synthesized without particle filter. Micrographs from left to right show respectively, the experimental high resolution micrograph, synthe- sized micrographs without particle filter which are statistically equivalent with high resolution micrograph at location C. High resolution micrograph at loca- tion C is used only for comparison and is not used in the synthesis process. . . 104 6.13 Micrograph at C synthesized with particle filter. From 6.13(a) to 6.13(c), exper- imental high resolution micrograph, experimental low resolution micrograph and high resolution micrograph synthesized with texture model and particle fil- ter at location C. High resolution micrograph at location C is used only for comparison and is not used in the synthesis process. . . . . . . . . . . . . . . . 105 xii Abstract This dissertation focuses on facilitating the analysis of probabilistic models for physical sys- tems. To that end, novel contributions are made to various aspects of the problem, namely, 1) development of efficient algorithms for solving stochastic system of equations governing the phys- ical problem, 2) stochastic basis adaptation methods to compute the solution in reduced stochas- tic dimensional space and 3) stochastic upscaling methods to find coarse-scale models from the fine-scale stochastic solutions. In particular, algorithms are developed for stochastic systems that are governed by partial differential equations (PDEs) with random coefficients. Polynomial chaos-based stochastic Galerkin and stochastic collocation methods are employed for solving these equations. Solvers and preconditioners based on Gauss-Seidel and Jacobi algorithms are explored for solving system of linear equations arising from stochastic Galerkin discretization of PDEs with random input data. Gauss-Seidel and Jacobi algorithms are formulated such that the existing software is leveraged in the computational effort. These algorithms are also used to de- velop preconditioners to Krylov iterative methods. These solvers and preconditioners are tested by solving a steady state diffusion equation and a steady state advection-diffusion equation. Upon discretization, the former PDE results in a symmetric positive definite matrix on left-hand-side, whereas the latter results in a non-symmetric positive definite matrix. The stochastic systems face significant computational challenge due the curse of dimensionality as the solution often lives in very high dimensional space. This challenge is addressed in the present work by recognizing the low dimensional structure of many quantities of interest (QoI) even in problems that have been embedded, via parameterization, in very high-dimensional settings. A new method for the char- acterization of subspaces associated with low-dimensional QoI is presented here. The probability xiii density function of these QoI is found to be concentrated around one-dimensional subspaces for which projection operators are developed. This approach builds on the properties of Gaussian Hilbert spaces and associated tensor product spaces. For many physical problems, the solution lives in multiple scales, and it is important to capture the physics at all scales. To address this issue, a stochastic upscaling methodology is developed in which the above developed algorithms and basis adaptation methods are used. In particular upscaling methodology is demonstrated by developing a coarse scale stochastic porous medium model that replaces a fine-scale which consists of flow past fixed solid inclusions. The inclusions have stochastic spatially varying thermal conductivities and generate heat that is transported by the fluid. The permeability and conductivity of the effective porous medium are constructed as statistically dependent stochastic processes that are both explicitly dependent on the fine scale random conductivity. Another contribution of this thesis is development of a probabilistic framework for synthesiz- ing high resolution micrographs from low resolution ones using a parametric texture model and a particle filter. Information contained in high resolution micrographs is relevant to the accurate prediction of microstructural behavior and the nucleation of instabilities. As these micrographs may be tedious and uneconomical to obtain over an extended spatial domain, A statistical ap- proach is proposed for interpolating fine details over a whole computational domain starting with a low resolution prior and high resolution micrographs available only at a few spatial locations. As a first step, a small set of high resolution micrographs are decomposed into a set of multi- scale and multi-orientation subbands using a complex wavelet transform. Parameters of a texture model are computed as the joint statistics of the decomposed subbands. The synthesis algorithm then generates random micrographs satisfying the parameters of the texture model by recursively xiv updating the gray level values of the pixels in the input micrograph. A density-based Monte Carlo filter is used at each step of the recursion to update the generated micrograph, using a low resolu- tion micrograph at that location as a measurement. The process is continued until the synthesized micrograph has the same statistics as those from the high resolution micrographs. The proposed method combines a texture synthesis procedure with a particle filter and produces good quality high resolution micrographs. xv Chapter 1 Introduction 1.1 Motivation Many physical problems can be described with coupled stochastic nonlinear partial differential equations. One such problem of interest is fluid flow through porous media coupled with heat equation. In this work a collection of two-dimensional solid discs is considered to be micro porous medium and the flow through pores at micro level can be modeled with Navier-Stokes equations coupled with heat equation. Due to inherent randomness in this physical system, it is modeled with uncertainty. Simulating such a system only at fine scale with Navier-Stokes equations can be computationally intractable. Hence, multi-scale modeling approach is employed where, coarse scale equations are modeled with Darcy-Brinkman equations. To avoid the com- plexity due to turbulence, the simulations are performed at a low Reynolds number, so that the flow is laminar. The coupling between the fluid flow and heat transfer is achieved through tem- perature dependent viscosity and fluid density. The multi-scale treatment of the problem makes the simulation numerically feasible but one has to be careful to maintain sufficient accuracy while making the multi-scale approximation. The main objective of this research is to come up with 1 stochastic upscaling approach so that the coarse scale porous media model can simulate the so- lution with required accuracy. Since the coarse scale models are much smaller compared to the fine-scale models, they are computationally less expensive to solve. Another problem considered in this work is probabilistic treatment of high resolution micro- graph synthesis from low resolution ones using a parametric texture model and a particle filter. High resolution microstructural information is required as a prerequisite for the multi-scale mod- eling and analysis of composite structures, and in particular to track the nucleation of instabilities in material behavior. As it remains expensive to experimentally obtain high resolution micro- graphs throughout a computational domain where coarse-scale behavior is being assessed, proce- dures for gleaning equivalent information from lower resolution images have demonstrated their value in computational materials science Ghosh et al. (2006), Valiveti and Ghosh (2007). 1.2 Problem statement In this work, stochastic upscaling methodology is developed to model the flow through porous media coupled with heat equations. At fine scale, fluid flow through an array of solid circular discs is modeled with Navier-Stokes equations Betchen et al. (2006), Chandesris et al. (2006), Nakayama et al. (2002) coupled with heat equation to model heat transfer between solid and fluid domain. Randomness in the the system is modeled by treating thermal conductivity of the solid region as a random filed. Spectral stochastic finite element methods are used to solve these stochastic par- tial differential equations. In these methods random coefficients and the solution are represented 2 Figure 1.1: Fine scale domain modeled with an array of solid discs and coarse scale domain modeled as porous media as series expansion of polynomial chaos basis, where, polynomial chaos basis are orthogonal polynomials in terms of input random variables. It is computationally intractable to solve Navier-Stokes equations with random parameters at the fine scale. To overcome this difficulty, an approximate porous media model is proposed at the coarse scale that will compute the solution of interest accurately. The equations governing the flow at coarse scale are modeled with Darcy-Brinkman equations Schmidt et al. (2010). A stochastic upscaling procedure is followed to compute the parameters of the coarse scale model. The parameters of the coarse scale model that are to be computed from the stochastic upscaling method are permeability and thermal conductivity random fields at the coarse scale. Figure 1.1 shows fine scale and coarse scale domain modeled in 2-dimensions. To synthesize high resolution micrographs, from low resolution ones, it is shown in chapter 6 that, a parametric texture model based in wavelet transformation and a particle filter can be used to accurately simulate the microstructural feature of the local domain. Starting with an image x 0 that represents the low-resolution micrograph at a location, our aim is to synthesize an image x at the same location, representing the high-resolution micrograph, that is constrained, in some 3 sense, by our knowledge of experimental high resolution data available at few locations and low resolution data available at the point of interest. In order to complete the statement of the problem, we must specify the nature of the constraints (i.e. the specific functionals or features that describe our knowledge of the high resolution data). In what follows, a texture model will be used to provide a context in which features of a micrograph are quantified. This model will be constructed from available high resolution data and presumed to be valid everywhere else. Specifically, the texture model will be identified with the joint statistics of the subbands in a particular wavelet decomposition, namely, a steerable pyramid decomposition. Associated with this model is a mapping f that maps x 0 into x such that x has the joint statistics corresponding to the model. In our work, we use a Bayesian framework where the experimental low resolution micrograph play the role of measurement used to update a prior model synthesized from a very small set of high resolution micrographs. We use a particle filter, namely a density-based Monte Carlo filter to implement the Bayesian formulation. 1.3 Stochastic upscaling Objective of stochastic upscaling in this work is to identify the coarse scale models that capture the quantities of interest with good accuracy. The parameters of the coarse scale models are computed using volume average method. Further, stochastic basis adaptation methods are developed to represent the coarse scale parameters in the reduced stochastic dimension. 4 1.4 Challenges There are many modeling, numerical and computational challenges involved in this problem. Significant modeling challenges in this problem are to choose the right governing equations and boundary conditions in both fine and coarse scale problems, to choose the spatial dependence of the coarse scale parameters and to come up with the appropriate distance measure between the quantities of interest at fine and coarse scales. Other modeling challenge is in representing the uncertainty at both the scales. The numerical challenges involved are instabilities in the finite element formulations of in- compressible flows. Main sources of these instabilities are spurious node to node oscillations in the velocity field in the advection dominated flow and spurious pressure modes in Navier-Stokes and Darcy-Brinkman equations due to equal order interpolation in pressure and velocity field in finite element discretizations. To address this issue, stabilization methods should be employed in the finite element discretization. In this work, stabilizations methods such as Streamline- upwind/Petrov-Galerkin (SUPG) and pressure-stabilizing/Petrov-Galerkin (PSPG) formulations Brooks and Hughes (1982), Shadid et al. (2006), Tezduyar (1992) are used. There are many computational challenges in solving this problem. At fine scale and coarse scale, the governing equations are coupled stochastic nonlinear partial differential equations. In the optimizations process, these partial differential equations should be solved efficiently and ac- curately. The stochastic partial differential equations at the fine scale are Navier-Stokes equations, which are very complex to solve even in the deterministic case. Spectral stochastic Galerkin pro- jections are used to discretize the stochastic partial differential equations. The linear algebraic equations arising from the discretization of these equations are very large in size and it is very 5 challenging to solve these equations efficiently. To address these issues, efficient solvers and preconditioners are proposed in chapter 3 based on Gauss-Seidel and Jacobi algorithms. To ex- plain Gauss-Seidel and Jacobi algorithms let Ax = b be a system of linear equations, where, A = [a ij ];b =fb i g andx =fx j g the Gauss-Seidel algorithms can be written as, x (k+1) j = 1 a ij 0 @ b i X j>i a ij x (k) j X j<i a ij x (k+1) j 1 A ; (1.1) and the Jacobi algorithm can be written as x (k+1) j = 1 a ij 0 @ b i X j6=i a ij x (k) j 1 A ; (1.2) These algorithms are adapted to solve the linear equations arising from spectral stochastic finite element discretization. The stochastic Galerkin system of equations for all the partial differential equations are formulated using Albany application software developed at Sandia National Labs. The solvers and preconditioners are incorporated as Stokhos package of Trilinos developed at Sandia National Labs. Another computational challenge arises due to curse of dimensionality in stochastic Galerkin projections. As the number of input random variables (stochastic dimension) increases, the num- ber of stochastic basis functions to represent the solution increases exponentially. Because of this curse of dimensionality, there is a limitation to go to larger stochastic dimensions. To address this issue by further reducing the size of the problem, dimension reduction methods based on basis adaptation discussed in chapter 4 will be used. In these methods, a new set of random variables in the reduced stochastic space is obtained b finding an isometry matrixA, the relates the new 6 set to the original set of random variables, = A. Here, =f 1 ; ; d g T are the original set of random variables and =f 1 ; ; l g T are the new set such that, l << d. The new set of random variables are adapted to specific quantities of interest (QoI) such that the solution of these QoI can be represented accurately with few random variables. Identifying the right set of random variables that are adapted to QoI is not unique. In chapter 4, various cases such as linear, quadratic basis adaptation is discussed. It is also shown that, if the QoI is scalar, it can be accurately represented in only on e random variable. Both, intrusive and non-intrusive methods can be employed for basis adaptation procedure. This manuscript is organized as follows. Chapter 2 provides a brief introduction to stochastic Galerkin methods. Chapter 3 describes the solvers and preconditiners based on Gauss-Seidel and jacobi algorithms to solve stochastic Galerkin system of equations efficiently. In chapter 4, ba- sis adaptation methods for dimension reduction is presented. Chapter 5 describes the stochastic upscaling methodology for a stochastic coupled nonlinear problem and shows the use of tools developed in chapters 3 and 4. Finally chapter 6 describes the probabilistic treatment of micro- graphs and their synthesis procedure from known statistical data. 7 Chapter 2 Stochastic Galerkin methods 2.1 Introduction Real life physical problems are often modeled as partial differential equations (PDEs) where the coefficients of the partial differential equations are treated as random parameters to represent uncertainty in the problem. Monte Carlo techniques are popular methods to solve these problems as they only require solutions to the PDE for a given set of realizations of the input random coefficients. More recently however, the stochastic finite element methods Babuˇ ska et al. (2004), Ghanem and Spanos (1991) have become a popular choice for solving these problems because of their advantages over Monte Carlo methods. These methods compute statistical properties of the solution more efficiently than Monte Carlo methods when the number of sources of uncertainty is not too large and the solution exhibits some degree of regularity with respect to these sources. 2.2 Stochastic partial differential equations LetD be an open subset ofR n and ( ; ;P ) be a complete probability space with sample space ,-algebra and probability measureP . We are interested in studying the following stochastic 8 partial differential equation: find a random field,u(x;!) : D !R such that the following holdsP -almost surely (P -a.s.): L(x;!;u) =f(x;!) in D ; (2.1) subject to the boundary condition B(x;!;u) =g(x;!) on@D ; (2.2) whereL is a differential operator andB is a boundary operator. 2.3 Input random field model Uncertainty in stochastic PDEs is often arises by treating coefficients in the differential operator, L(x;!;u), as a random fields. For computational purposes, each stochastic coefficienta(x;!) must be discretized in both spatial and stochastic domains. To this end, it is often approximated with a truncated series expansion that separates the spatial variablex from the stochastic variable ! resulting in a representation by a finite number of random variables. For this representation, second order information such as the covariance function of the random field is required. In the present problem, two cases of random field models are considered. In the first case, the random field is assumed to be uniformly distributed and is approximated through a truncated Karhunen-Lo` eve expansion. In the second case, the random field is assumed to have a log-normal distribution, that is a(x;!) = exp (g(x;!)) where g(x;!) is a Gaussian random field, and is approximated by a truncated polynomial chaos expansion. 9 2.3.1 Karhunen-Lo` eve expansion Let C(x 1 ;x 2 ) = E[a(x 1 ;!)a(x 2 ;!)] be the covariance function of the random field a(x;!), where E[] denotes mathematical expectation. Then a(x;!) can be approximated through its truncated Karhunen-Lo` eve (K-L) expansion Ghanem and Spanos (1991) given by a(x;!) ~ a(x;(!)) =a 0 (x) + M X i=1 p i a i (x) i (!); (2.3) where a 0 (x) is the mean of the random field a(x;!) andf( i ;a i (x))g i1 are solutions of the integral eigenvalue problem Z D C(x 1 ;x 2 )a i (x 2 )dx 2 = i a i (x 1 ): (2.4) The eigenvalues i are positive and non-increasing, and the eigenfunctionsa i (x) are orthonormal, that is, Z D a i (x)a j (x) = ij ; (2.5) where ij is the Kronecker delta. In Eq. 2.3,f i g M i=1 are uncorrelated random variables with zero mean. As a first test-case, the diffusion coefficient a(x;!) is modeled as a random field with following exponential covariance function C(x 1 ;x 2 ) = 2 exp(kx 1 x 2 k 1 =L) (2.6) such that the random variables i (!) in the K-L expansion are uniformly distributed. We further assume that the random variables are independent. 10 2.3.2 Polynomial chaos expansion The K-L expansion approximates a random field by a linear combination of a finite set of random variables. To maintain positivity of the random field, such a representation is only appropriate if the random variables are bounded Ullmann (2010). For unbounded random variables (e.g., log- normal) a nonlinear polynomial chaos representation is more appropriate. The polynomial chaos expansion Ghanem and Spanos (1991), Weiner (1963) is used to approximate a random field in terms of multi-variate orthogonal polynomials. Let = ( 1 ; ; M ) T be the random variables from a truncated K-L expansion of a given random fieldg(x;!), that is g(x;!) ~ g(x;(!)) =g 0 (x) + M X i=1 p i g i (x) i (!): (2.7) Assumea(x;!) is then given by a nonlinear transformation ofg(x;!). Thena(x;!) can be repre- sented through nonlinear functionals of the random variables i (!). It has been shown in Ghanem and Spanos (1991), Weiner (1963) that this functional dependence can be expanded in terms of multi-dimensional orthogonal polynomials, called polynomial chaos, as a(x;!) = ^ a 0 (x) + 1 X i 1 =1 ^ a i 1 (x) 1 ( i 1 (!)) + 1 X i 1 =1 i 1 X i 2 =1 ^ a i 1 i 2 2 ( i 1 (!); i 2 (!)) + (2.8) where n ( i 1 ; ; in ) is the multi-dimensional polynomial chaos of ordern in random variables ( i 1 ; ; in ). A one-to-one mapping of polynomialsf i g to a set of polynomials with ordered indicesf i ()g can be introduced Ghanem and Spanos (1991). After substitutingf i g in Eq. 2.8 11 and truncating the series to finite number of terms N , the random field a(x;!) can thus be approximated as a(x;!) ~ a(x;(!)) =a 0 (x) + N X i=1 a i (x) i (): (2.9) The polynomialsf i ()g are orthogonal with respect to the inner product defined by expectation in the stochastic space, h i (); j ()i Z i ((!)) j ((!))dP (!) = ij : (2.10) As a second test-case, the diffusion coefficienta(x;!), is modeled as a log-normal random field Ghanem (1999) wherea(x;!) = exp[g(x;!)] andg(x;!) is a Gaussian random field with exponential covariance (2.6). Here,g(x;!) is approximated with a truncated K-L expansion (2.7). In this case the random variables i are standard normal random variables and thus are inde- pendent and a(x;!) can be approximated with truncated polynomial chaos expansion (2.9). It also can be shown that the polynomialsf i g are tensor products of one-dimensional Hermite polynomials. For a given total polynomial orderp, the total number of polynomialsf i ()g is N + 1 = (M+p)! M!p! . 2.4 Stochastic Galerkin method LetH 1 0 (D) be the subspace of the Sobolev spaceH 1 (D) that vanishes on the boundary@D and is equipped with the normkuk H 1 0 (D) = [ R D jruj 2 dx] 1 2 . Problem (2.1) can then be written in the 12 following equivalent variational form Ghanem and Doostan (2006): findu2 H 1 0 (D) L 2 ( ) such that b(u;v) =l(v); 8v2H 1 0 (D) L 2 ( ); (2.11) whereb(u;v) is a continuous and coercive bilinear form andl(v) is a continuous bounded lin- ear functional. In the stochastic Galerkin method, we seek the solution of the variational prob- lem (2.11) in a tensor product spaceX h Y p , where,X h H 1 0 (D) is finite dimensional space of continuous polynomials corresponding to the spatial discretization ofD andY p L 2 ( ) is the space of random variables spanned by polynomial chaos Ghanem and Spanos (1991) of order up top. Then the finite dimensional approximationu X h Yp (x;!) of the exact solutionu(x;!) on the tensor product spaceX h Y p is given as the solution to b(u X h Yp ;v) =l(v) 8v2X h Y p : (2.12) In equation (2.12) the input random fielda(x;!) in the bilinear formb(u X h Yp ;v) can be ap- proximated using either a K-L expansion or a polynomial chaos expansion depending on a choice of model for the random field. The finite dimensional approximationu X h Yp (x;!) is represented with truncated polynomial chaos expansion, where the multidimensional polynomial chaos are orthogonal with respect to the probability measure of the underlying random variables. The re- sulting set of coupled PDEs are then discretized using standard techniques such as the finite ele- ment or finite difference methods. In the present problem, the set of coupled PDEs are discretized using the finite element method and the resulting system of linear equations can be written as 13 N X j=0 ^ P X i=0 c ijk K i u j =f k ; k = 0; ;N ; (2.13) wheref k =Eff(x;) k g,c ijk =Ef i j k g and ^ P =M whena(x;!) is approximated by a truncated K-L expansion, orc ijk = Ef i j k g and ^ P = ^ N whena(x;!) is approximated by a polynomial chaos expansion. HerefK i 2R NxNx g ^ P i=0 are the polynomial chaos coefficients of the stiffness matrix (section (3.4) of Powell and Elman (2009)) andfu j 2 R Nx g N j=0 are the polynomial chaos coefficients of the discrete solution vector u j = [u 0j ;:::;u Nxj ] T ; j = 0;:::;N : (2.14) Equation (2.13) can be written in the form of a global stochastic stiffness matrix of size ((N + 1)N x ) by ((N + 1)N x ) as 2 6 6 6 6 6 6 6 6 6 6 4 K 0;0 K 0;1 K 0;N K 1;0 K 1;1 K 1;N . . . . . . . . . . . . K N ;0 K N ;1 K N ;N 3 7 7 7 7 7 7 7 7 7 7 5 8 > > > > > > > > > > < > > > > > > > > > > : u 1 u 2 . . . u N 9 > > > > > > > > > > = > > > > > > > > > > ; = 8 > > > > > > > > > > < > > > > > > > > > > : f 1 f 2 . . . f N 9 > > > > > > > > > > = > > > > > > > > > > ; (2.15) whereK j;k = P ^ P i=0 c ijk K i . We will denote this system as K u = f. In practice it is prohibitive to assemble and store the global stochastic stiffness matrix in this form, rather each block of the stochastic stiffness matrix can be computed from thefK i g when needed. The stochastic stiff- ness matrix is block sparse, which means that some of the off-diagonal blocks are zero matrices, because of the fact that the c ijk defined above vanishes for certain combination of i;j and k. 14 It is also interesting to note that the diagonal blocks are dominant over the off-diagonal blocks because, the diagonal blocks have the contribution from mean stiffness matrixK 0 , whereas the off-diagonal blocks do not have the contribution from mean stiffness matrix. These properties of block sparsity and diagonal dominance are exploited to develop solvers and preconditioners based on Jacobi and Gauss-Seidel algorithms. 15 Chapter 3 Iterative solvers and preconditioners 3.1 Introduction The linear algebraic equations arising from the discretization of these equations are very large in size and it is very challenging to solve these equations efficiently. In this work solvers and preconditioners based on Jacobi and Gauss-Seidel algorithms are developed to solve the system of equations arising from the stochastic finite element methods. 3.2 Solution methods for stochastic Galerkin systems In this section, various solver techniques and preconditioning methods for solving the linear algebraic equations arising from stochastic Galerkin discretizations (2.13) are described. The solver methods discussed are: relaxation methods, namely, a Jacobi mean method and a Gauss- Seidel mean method, and Krylov-based iterative methods Saad (1996). Also various stochastic preconditioners used to accelerate convergence of the Krylov methods are discussed, including mean-based Powell and Elman (2009), Gauss-Seidel mean, approximate Gauss-Seidel mean, ap- proximate Jacobi mean and Kronecker product Ullmann (2010) preconditioners. The relaxation 16 schemes can be viewed as fixed point iterations on a preconditioned system Saad (1996). In Jacobi and Gauss-Seidel methods, mean splitting is used rather than traditional diagonal block splitting as it allows use of the same mean matrixK 0 for all inner deterministic solves (and thus reuse of the preconditionerP 0 K 0 ). 3.2.1 Jacobi mean algorithm In this method, systems of equations of size equal to that of the deterministic system are solved iteratively by updating the right-hand-side to obtain the solution to the stochastic Galerkin system of equations (2.13): c kk0 K 0 u new k =f k N X j=0 ^ P X i=1 c ijk K i u old j ; k = 0; ;N : (3.1) The above system of equations are solved fork = 0; ;N using any solution technique appro- priate for the mean matrixK 0 . Thus existing legacy software can be used with minimal modifi- cation to solve the stochastic Galerkin system. In this work, Krylov-based iterative methods with appropriate preconditioners will be used. One cycle of solves fromk = 0; ;N is considered one Jacobi outer iteration, and after each outer iteration, the right-hand-side in equation (3.1) is updated replacingfu old j g with the new solutionfu new j g. These outer iterations are continued until the required convergence tolerance is achieved. The Jacobi mean algorithm is shown in Algorithm 1. 17 Algorithm 1 Jacobi mean algorithm 1. Choose initial guess u 0 and compute residual r = K u 0 f 2. Iteration count,itr = 0 3. while k rk 2 k fk 2 > tol do 4. fork = 0:::N do 5. Solvec kk0 K 0 u (itr+1) k =f k P N j=0 P ^ P i=1 c ijk K i u (itr) j 6. end for 7. itr =itr + 1 8. r = K u itr f 9. end while. Note that for a given outer iteration, all of the right-hand-sides fork = 0; ;N are available simultaneously, and thus their solution can be efficiently parallelized. Moreover block algorithms optimized for multiple right-hand-sides may be used to further increase performance. Finally this approach does not require a large amount of memory to compute the solution. The disadvantage of the method is that it may not converge or may converge very slowly when the diagonal blocks of the stochastic stiffness matrix are less dominant over off-diagonal blocks. 3.2.2 Gauss-Seidel mean iterative method The Gauss-Seidel method considered is similar to the Jacobi method, except the right-hand-side in equation (3.1) is updated after each deterministic solve with the newly computedu new k . Sym- bolically this is written c kk0 K 0 u new k =f k k1 X j=0 ^ P X i=1 c ijk K i u new j N X j=k ^ P X i=1 c ijk K i u old j ; k = 0; ;N : (3.2) 18 As before, one cycle of solves fromk = 0; ;N is considered one outer iteration of the Gauss- Seidel method, and these outer iterations are repeated until the required convergence tolerance is achieved. Note however that computing the updates as shown here would result in a large number of duplicated matrix-vector products K i u j for each outer iteration. Instead, after each u new k is computed by solving the mean linear system, we first computey =K i u new k for alli in whichc ijk is nonzero for anyj. Then for each correspondingj we updatef j f j c ijk y. This allows all of the right-hand-sides to be updated as required using the fewest number of matrix-vector products and without resorting to storing intermediate products. The complete Gauss-Seidel algorithm is shown in Algorithm 2. 19 Algorithm 2 Gauss-Seidel mean algorithm 1. Choose initial guess u 0 and compute residual r = f K u 0 2. Iteration count,itr = 0 3. Initialize z = r 4. while k rk 2 k fk 2 > tol do 5. r = f 6. fork = 0:::N do 7. Solvec kk0 K 0 u k =z k 8. z k =f k 9. fori = 1;:::; ^ P 10. y =K i u k 11. forj = 1;:::;N 12. ifc ijk 6= 0 then 13. z j =z j c ijk y 14. r j =r j c ijk y 15. endif 16. end for 17. end for 18. r k =r k c kk0 K 0 u k 19. end for 20. itr =itr + 1 21. end while. 20 Often this method converges in fewer iterations than the Jacobi method, at the expense of no longer having all of the right-hand-sides available simultaneously. Unlike diagonal block splitting methods defined in Rosseel and Vandewalle (2010), mean block splitting is used in both Jacobi and Gauss-Seidel algorithms and hence the left-hand-side matrix is the mean matrix for all inner deterministic problems and only the right-hand-side changes. In such cases recycled Krylov basis methods could be explored to increase performance. 3.2.3 Krylov based iterative methods with matrix-free operations Krylov based iterative methods Saad (1996) such as the conjugate gradient (CG) method and the generalized minimal residual (GMRES) method can be used to solve the stochastic Galerkin sys- tem (2.13) in which matrix vector products v = K u are computed using “matrix-free” operations: v k = N X j=0 ^ P X i=0 c ijk K i u j ; k = 0; ;N : (3.3) If the matrix vector products are computed from Eq. 3.3, it is not required to assemble the full stochastic Galerkin stiffness matrix, drastically decreasing memory requirements. However if a large number of iterations of a Krylov method such as GMRES are required, allocation of the Krylov basis may still require a very large amount of memory. Thus good preconditioning strategies for the stochastic Galerkin system are required, several of which will be discussed below. 21 Mean-based preconditioner The mean-based preconditioner Powell and Elman (2009) is given by P = diagfP 0 ; ; P 0 g whereP 0 K 0 is a preconditioner for the mean. The mean-based preconditioner is very efficient to compute and apply, since it only must be generated once from a matrix that is of the size of the deterministic system. However it doesn’t incorporate any higher-order stochastic information, thus its performance degrades as the stochastic dimension, polynomial order, or random field variance increases Ullmann (2010). Gauss-Seidel preconditioner One or more outer iterations of the Gauss-Seidel mean algorithm can be used as a preconditioner to the Krylov based iterative methods. An advantage of this method is that the cost of applying the preconditioner can be controlled by adjusting the tolerance of the inner deterministic solves and number of outer iterations. Decreasing this tolerance and increasing the number of outer it- erations will reduce the number of iterations in the Krylov method, but make the preconditioner more expensive to apply, and thus these must be balanced to minimize overall computational cost. Generally we have found the cost of the preconditioner to be dominated by solving the mean sys- tems, and thus the performance was improved by loosening the outer solver tolerance or limiting the number of outer iterations. In the results presented below we limited the preconditioner is limited to only one Gauss-Seidel iteration. 22 Approximate Gauss-Seidel preconditioner The process of increasing the inner solver tolerance can be taken to its extreme of replacing the inner mean solves by application of the mean preconditioner. As with the Gauss-Seidel precondi- tioner above, we found experimentally that this approach worked best with only one Gauss-Seidel iteration, and adding additional iterations did not improve the quality of the preconditioner. We also found that the cost of the preconditioner was reduced dramatically if only the first-order terms in the expansion for the stiffness matrix were used in the preconditioner and using higher-order terms did not improve performance. It is referred as the approximate Gauss-Seidel preconditioner here. Approximate Jacobi preconditioner Similar to the approximate Gauss-Seidel preconditioner, Jacobi iterations can be used using a preconditioner in place of the mean stiffness matrix. In this case we used two outer Jacobi it- erations, since the first iteration is equivalent to mean-based preconditioning (i.e., the additional terms on the right-hand-side of equation (3.1) are zero). Increasing the number of outer iterations did not improve the efficiency of the overall solver. It is referred as the approximate Jacobi pre- conditioner. Both approximate Gauss-Seidel and approximate Jacobi preconditioners are found to be very effective in reducing the number of Krylov iterations and are also not very expensive to apply. 23 Kronecker product preconditioner The Kronecker product preconditioner Ullmann (2010) is defined asP 1 =G K 0 , whereK 0 is mean stiffness matrix andG is G = ^ P X i=0 tr(K T i K 0 ) tr(K T 0 K 0 ) G i (3.4) where,G i (j;k) = c ijk . Unlike, the mean-based preconditioner, the Kronecker product precon- ditioner incorporates higher order stochastic information and hence the Krylov based algorithm converges in fewer iterations. However the disadvantage is that the cost of constructing the Kro- necker product preconditioner is larger than that of mean-based preconditioner, and it is also more expensive to apply. Over all, solution time is found to be less than that with mean-based precon- ditioner as reported in Ullmann (2010). However, it was found in the numerical experiments that the approximate Gauss-Seidel and approximate Jacobi preconditioner performed better than the Kronecker product preconditioner. 3.3 Problem Statement In this work a stochastic steady-state elliptic diffusion equation with Dirichlet boundary condi- tions Elman et al. (2011) is used as test problem for various stochastic PDE solution methods. Assumea(x;!) :D !R is a random field that is bounded and strictly positive, that is, 0<a l a(x;!)a u <1 a:e: in D : (3.5) 24 For the steady-state diffusion equation, we wish to compute a random fieldu(x;!) :D !R, such that the following holdsP -almost surely (P -a.s.): r:(a(x;!)ru(x;!)) =f(x;!) in D ; (3.6) u(x;!) =u 0 on@D : (3.7) Problem (3.6) can then be written in the following equivalent variational form Ghanem and Doostan (2006): findu2H 1 0 (D) L 2 ( ) such that b(u;v) =l(v); 8v2H 1 0 (D) L 2 ( ); (3.8) whereb(u;v) is the continuous and coercive (from assumption (3.5)) bilinear form given by b(u;v) =E Z D arurvdx ; 8u;v2H 1 0 (D) L 2 ( ): (3.9) andl(v) is the continuous bounded linear functional given by l(v) =E Z D fvdx ; 8v2H 1 0 (D) L 2 ( ): (3.10) From the Lax-Milgram lemma, equation (3.8) has unique a solution inH 1 0 (D) L 2 ( ). 25 3.4 Sparse grid collocation method In the collocation method, the solution to the PDE is sampled at a pre-selected set of points called collocation points, = ( (1) ; ; (N) ). The stochastic solution is constructed by interpolating at these collocation points, u(x;) N X k=0 u k (x)L k (); (3.11) wherefL k ()g are Lagrange interpolatory polynomials defined by (L k ( l ) = kl ) andu k is the solution of following deterministic diffusion equation, r:(a(x; (k) )ru k (x) =f(x; (k) ) in D; (3.12) u k (x) =u k0 on@D: (3.13) The collocation points can be chosen as tensor products of 1-D Gaussian quadrature points and the interpolating polynomials as tensor products of 1-D Lagrange interpolating functions. However the number of collocation points then grows exponentially with the number of random variables. An alternative method is to use Smolyak sparse grid quadrature (Nobile et al. (2008), Smolyak (1963)) where collocation points that do not increase asymptotic accuracy are removed from the tensor product grid. This results in many fewer collocation points but still more than the number of stochastic degrees of freedom in the stochastic Galerkin method. This method is fully non-intrusive and is easy to implement with existing legacy software (once the sparse grid is generated) Adams et al. (2010). 26 3.5 Numerical illustration To compare the performance of different solvers and preconditioners discussed above, a 2-D stochastic diffusion equation and a 2-D stochastic advection-diffusion equation presented in sec- tion 3.3 are solved using the stochastic Galerkin method described in section 2.4. In the above two problems, the diffusion coefficient is modeled as both a random field discretized using a truncated K-L expansion with uniform random variables (section 2.3.1) and a log-normal random field discretized using a truncated polynomial chaos expansion (section 2.3.2). In the first case, the orthogonal polynomials used in the stochastic Galerkin method are tensor products of 1-D Legendre polynomials and in the second case tensor products of Hermite polynomials are used. For simplicity a constant unit forcef(x;!) = 1 is used as the right-hand-side in equation (3.6). The spatial dimensions are discretized using standard finite element mesh with linear quadrilat- eral elements. In advection-diffusion equation, the parameter, ~ w = [1; 1] T . The corresponding stochastic Galerkin linear system is constructed using the Stokhos and Epetra packages in Trilinos. For the Jacobi solver, Gauss-Seidel solver, Gauss-Seidel preconditioner, and stochastic colloca- tion method, the linear systems are solved via multi-grid preconditioned GMRES provided by the AztecOO and ML packages in Trilinos. To use CG, a symmetric version of the Gauss-Seidel pre- conditioner is required which comprises of one forward and one backward Gauss-Seidel iteration, and hence it is twice as expensive as that of one forward Gauss-Seidel iteration. It was observed from numerical tests that the symmetric Gauss-Seidel improved the number of Krylov iterations slightly, though its cost is twice that of the non-symmetric Gauss-Seidel preconditioner. Hence, GMRES algorithm is employed for both diffusion and advection-diffusion problems instead of CG. 27 For the numerical comparisons, the domainD = [0:5; 0:5] [0:5; 0:5] is discretized into a 64 64 grid resulting in a total number of nodes,N x = 4096. In figure (3.1) the solution time for the stochastic Galerkin method, scaled by the deterministic solution time at the mean of the random field, is compared for different mesh sizes (32 32, 64 64, 96 96 and 128 128) for the diffusion problem demonstrating that the solution time does not depend strongly on the mesh size (as is to be expected for the multigrid preconditioner). The scaled solution time for these solvers and preconditioning techniques as a function of the standard deviation of the input random field, stochastic dimension, and polynomial order are then tabulated in Tables 3.1-3.12. In the tables, the number of Krylov iterations for the aforementioned preconditioners and iterations for Gauss-Seidel and Jacobi solvers are provided in parentheses. In the tables, MB, AGS, AJ, GS and KP are the mean-based, approximate Gauss-Seidel, approximate Jacobi, Gauss-Seidel and Kronecker-product preconditioners respectively. GS in the solution methods refers to the Gauss- Seidel mean algorithm (Algorithm-2). “Jacobi” refers to the Jacobi mean algorithm (Algorithm- 1). The solution tolerance for all of the stochastic Galerkin solvers is 1e 12 . For the Gauss-Seidel and Jacobi solvers, the inner solver tolerance is 3e 13 . All the computations are performed using a single core of an 8 core, Intel Xeon machine with 2.66 GHz and 16GB Memory. Figures 3.2 and 3.3 show the plots of relative residual error vs iteration count for the stochastic Galerkin system with stochastic dimension 4 and polynomial order 4 and standard deviation 0.1. It can be observed that the matrix-free Krylov solver with the Gauss-Seidel preconditioner takes the least number of iterations in case of uniform random field and Gauss-Seidel solver in case of log-normal random field, whereas the Jacobi solver takes highest number of iterations in both cases for a given tolerance. However in terms of solution time, the matrix-free Krylov solver with the approximate Gauss-Seidel preconditioner is the most efficient compared to all other stochastic 28 Galerkin solvers. Comparison in terms of iteration count alone is misleading to evaluate the computational cost because, in each iteration the cost of preconditioner as observed in the case of Gauss-Seidel preconditioner could be very high resulting in higher computational cost even with small number of iterations. Hence we also compare the solution time for all preconditioners and solvers. In the tables, “Div” means diverged. In the case of the uniform random field with small variance ( = 0:1), it can be observed from Tables 3.1 and 3.2 that more intrusive Krylov-based stochastic Galerkin solvers are more efficient than less intrusive Gauss-Seidel and Jacobi solvers. Moreover the approximate Gauss-Seidel and Jacobi preconditioners are a significant improvement over the traditional mean-based approach. However as the variance of the random field increases, we can see from Table 3.3 that the Gauss-Seidel and Jacobi solvers suffer considerably, whereas the the Krylov-based approaches (excluding the Gauss-Seidel preconditioner) still perform quite well. This is not unexpected, as the operator becomes more indefinite as the variance increases. In the case of the advection-diffusion equation, the stiffness matrix is non-symmetric. Tables 3.4-3.6 show the scaled solution time for the advection-diffusion equation with the diffusion coefficient modeled using a Karhunen-Lo` eve with uniform random variables. We can observe that the Krylov solver with the approximate Gauss-Seidel preconditioner performs better than the rest of the Galerkin solvers. In the case of the random field with uniform random vari- ables, the solution time for the diffusion equation is smaller than that for the advection-diffusion equation. This observation could be due to the fact that the multi-grid preconditioner for solving the deterministic diffusion equation is more efficient than that for solving the advection-diffusion equation. Since the efficiency of the preconditioners and solvers based on Gauss-Seidel and Jacobi algorithms depends heavily on the efficiency of the deterministic solver, a good preconditioner for 29 the mean stiffness matrix will significantly improve the above mentioned stochastic Galerkin pre- conditioners and solvers. In the case of the log-normal random field, we can see from Tables 3.7 and 3.11 that Gauss- Seidel and Jacobi solvers have not performed well in terms of solution time. The Gauss-Seidel solver is comparable to the Krylov solver with the AGS preconditioner only in case of log-normal random field and at higher polynomial chaos order. For higher variance of the random field, we see from Table 3.9 that the Jacobi solver diverges. This problem might be addressed to some extent by using the true diagonal matrixK k;k = P M i=0 c ikk K i from global stochastic stiffness matrix as the left-hand-side in the Jacobi solver and preconditioner instead of the mean matrixK 0 . In the case of the log-normal random field, both the diffusion and advection-diffusion equations have similar trends and the solution times are comparable. From the Krylov iteration count provided in the parentheses of the tables, we can observe that the preconditiners and solvers based on Gauss- Seidel and Jacobi are robust with respect to stochastic dimension and polynomial order, however they are not robust with respect to the variance of the input random field. (a) Stochastic dimension = 3 (b) Stochastic dimension = 4 Figure 3.1: Scaled Galerkin solution time vs solution error for varying polynomial chaos order and spatial mesh for the diffusion problem with uniform random field 30 (a) difusion coefficient modeled as uniform ran- dom field (b) difusion coefficient modeled as lognormal random field Figure 3.2: Relative residual norm vs iteration countto solve stochastic diffusion equation (a) difusion coefficient modeled with uniform random variables (b) difusion coefficient modeled as lognormal random field Figure 3.3: Relative residual norm vs iteration count to solve stochastic advection-diffusion equa- tion 31 Table 3.1: Scaled solution time (# of iterations) vs stochastic dimension for random field with uniform random variables, diffusion equation, PC order = 4 and=0.1 Stoch. Preconditioners for GMRES GS, Jacobi Solvers dim MB AGS AJ GS KP GS Jacobi 1 10 (21) 6.5 (13) 12 (13) 46 (9) 7 (14) 67 (15) 121 (29) 2 39 (25) 24 (15) 43 (15) 168 (11) 28 (18) 256 (19) 490 (39) 3 108 (28) 62 (16) 109 (16) 429 (12) 86 (22) 723 (23) 1355 (46) 4 235 (29) 138 (17) 238 (17) 937 (13) 188 (23) 1587 (25) 2973 (50) 5 482 (30) 267 (17) 454 (17) 1706 (13) 408 (25) 2980 (26) 5725 (53) 6 890 (32) 490 (18) 822 (18) 3093 (14) 776 (27) 5639 (29) 10988 (59) 7 1456 (33) 787 (18) 1308 (18) 4902 (14) 1241 (27) 9361 (30) 19666 (61) 32 Table 3.2: Scaled solution time (# of iterations) vs order of polynomial chaos for random field with uniform random variables, diffusion equation, Stoch. dim=4, = 0:1 PC Preconditioners for GMRES GS, Jacobi Solvers order MB AGS AJ GS KP GS Jacobi 2 34 (22) 22 (14) 40 (14) 140 (9) 28 (18) 215 (16) 408 (32) 3 99 (26) 62 (16) 103 (15) 397 (11) 81 (21) 660 (21) 1216 (41) 4 235 (29) 138 (17) 238 (17) 937 (13) 188 (23) 1587 (25) 2973 (50) 5 523 (32) 285 (18) 480 (18) 1850 (14) 411 (25) 3344 (29) 6366 (59) 6 948 (35) 516 (19) 863 (19) 3265 (15) 774 (27) 6464 (32) 12327 (67) 7 1582 (37) 821 (19) 1439 (20) 5536 (16) 1338 (29) 11354 (35) 24511 (73) Table 3.3: Scaled solution time (# of iterations) vs standard deviation () of input random field, for random field with uniform random variables, diffusion equation, Stoch dim = 4 and PC order = 4 Preconditioners for GMRES GS, Jacobi Solvers MB AGS AJ GS KP GS Jacobi 0.10 235 (29) 138 (17) 238 (17) 937 (13) 188 (23) 1587 (25) 2973 (50) 0.11 272 (33) 147 (18) 255 (18) 1005 (14) 217 (26) 1898 (30) 3666 (61) 0.12 298 (37) 164 (20) 283 (20) 1139 (16) 245 (29) 2374 (37) 4585 (75) 0.13 330 (41) 182 (22) 316 (22) 1340 (19) 283 (33) 3057 (48) 5998 (97) 0.14 384 (48) 211 (25) 361 (25) 1544 (22) 318 (38) 4152 (65) 8435 (133) 0.15 487 (60) 272 (31) 458 (31) 1956 (28) 382 (46) 6567 (101) 13525 (206) 33 Table 3.4: Scaled solution time (# of iterations) vs stochastic dimension for random field with uniform random variables, advection-diffusion equation, PC order = 4, and=0.1 Stoch. Preconditioners for GMRES GS, Jacobi Solvers dim MB AGS AJ GS KP GS Jacobi 1 12 (24) 7 (15) 14 (15) 65 (12) 9 (18) 64 (14) 114 (27) 2 47 (30) 27 (17) 48 (17) 240 (15) 37 (23) 256 (19) 489 (39) 3 133 (34) 74 (19) 129 (19) 631 (17) 106 (27) 727 (23) 1353 (46) 4 291 (36) 165 (20) 285 (20) 1343 (18) 246 (29) 1579 (25) 3088 (52) 5 598 (38) 339 (21) 562 (21) 2571 (19) 532 (31) 3098 (27) 5941 (55) 6 1068 (40) 611 (22) 1012 (22) 4994 (21) 989 (34) 5808 (30) 11567 (62) 7 1740 (41) 976 (22) 1690 (22) 7589 (21) 1619 (35) 9951 (32) 21199 (64) 34 Table 3.5: Scaled solution time (# of iterations) vs order of polynomial chaos for random field with uniform random variables, advection-diffusion equation, Stoch. dim=4, = 0:1 PC Preconditioners for GMRES GS, Jacobi Solvers order MB AGS AJ GS KP GS Jacobi 2 41 (26) 26 (16) 45 (16) 213 (13) 35 (22) 233 (17) 429 (33) 3 125 (32) 70 (18) 124 (18) 604 (16) 102 (26) 663 (21) 1278 (43) 4 291 (36) 165 (20) 285 (20) 1343 (18) 246 (29) 1579 (25) 3088 (52) 5 626 (40) 336 (21) 591 (22) 2654 (20) 543 (32) 3302 (29) 6436 (60) 6 1163 (44) 637 (23) 1104 (24) 5082 (22) 1008 (35) 6311 (33) 12558 (68) 7 1981 (47) 1062 (24) 1825 (25) 8443 (23) 1713 (38) 11143 (36) 25117 (74) Table 3.6: Scaled solution time (# of iterations) vs standard deviation () for random field with uniform random variables, advection-diffusion equation, Stoch dim = 4, PC order = 4 Preconditioners for GMRES GS, Jacobi Solvers MB AGS AJ GS KP GS Jacobi 0.10 291 (36) 165 (20) 285 (20) 1343 (18) 246 (29) 1579 (25) 3088 (52) 0.11 328 (41) 182 (22) 326 (23) 1556 (21) 284 (33) 1960 (31) 3711 (62) 0.12 389 (49) 218 (26) 371 (26) 1921 (25) 325 (39) 2404 (38) 4624 (76) 0.13 469 (58) 257 (30) 461 (31) 2430 (30) 380 (46) 3039 (48) 5975 (97) 0.14 589 (73) 324 (38) 565 (39) 2968 (38) 478 (57) 4123 (65) 8153 (131) 0.15 810 (101) 437 (52) 742 (52) 3949 (52) 650 (78) 6110 (96) 12514 (196) 35 Table 3.7: Scaled solution time (# of iterations) vs stochastic dimension for log-normal random field, diffusion equation, PC order = 4, and=0.1 Stoch. Preconditioners for GMRES GS, Jacobi Solvers dim MB AGS AJ GS KP GS Jacobi 1 9 (16) 8 (13) 12 (12) 48 (9) 7 (12) 35.87 (8) 86 (20) 2 33 (17) 26 (13) 40 (12) 149 (9) 28 (14) 108 (8) 255 (20) 3 90 (17) 72 (13) 102 (12) 363 (9) 81 (15) 254 (8) 629 (21) 4 231 (17) 183 (13) 244 (12) 782 (9) 217 (15) 534 (8) 1327 (21) 5 594 (17) 473 (13) 577 (12) 1603 (9) 580 (15) 922 (8) 2292 (21) 6 1348 (17) 1060 (13) 1220 (12) 3060 (9) 1387 (15) 1613 (8) 4236 (21) Table 3.8: Scaled solution time (# of iterations) vs order of polynomial chaos for log-normal random field, diffusion equation, Stoch. dim=4, = 0:1 PC Preconditioners for GMRES GS, Jacobi Solvers order MB AGS AJ GS KP GS Jacobi 2 25 (15) 23 (13) 36 (12) 116 (7) 24 (14) 94 (7) 202 (16) 3 77 (16) 65 (13) 97 (12) 355 (9) 70 (14) 253 (8) 566 (19) 4 231 (17) 183 (13) 244 (12) 782 (9) 217 (15) 534 (8) 1327 (21) 5 738 (18) 593 (14) 644 (12) 1888 (10) 664 (15) 918 (8) 2642 (24) 6 2181 (19) 1645 (14) 1660 (12) 4111 (10) 1963 (15) 1606 (8) 5046 (26) 36 Table 3.9: Scaled solution time (# of iterations) vs standard deviation () for log-normal random field, diffusion equation, Stoch dim = 4 and PC order = 4 Preconditioners for GMRES GS, Jacobi Solvers MB AGS AJ GS KP GS Jacobi 0.10 231 (17) 183 (13) 244 (12) 782 (9) 217 (15) 534 (8) 1327 (21) 0.15 275 (20) 217 (15) 267 (13) 942 (11) 243 (17) 638 (10) 2030 (33) 0.20 332 (24) 247 (17) 285 (14) 1184 (14) 273 (19) 764 (12) 3897 (62) 0.25 377 (27) 286 (20) 288 (14) 1340 (16) 302 (21) 957 (15) 18672 (268) 0.30 456 (32) 328 (23) 346 (17) 1583 (19) 333 (23) 1150 (18) Div 0.35 518 (37) 372 (26) 591 (29) 1827 (22) 359 (25) 1406 (22) Div Table 3.10: Sacled solution time (# of iterations) vs stochastic dimension for log-normal random field, advection-diffusion equation, PC order = 4 and=0.1 Stoch. Preconditioners for GMRES GS, Jacobi Solvers dim MB AGS AJ GS KP GS Jacobi 1 9 (16) 8 (13) 12 (12) 48 (9) 12 (12) 36 (8) 81 (19) 2 33 (17) 26 (13) 40 (12) 148 (9) 40 (12) 108 (8) 255 (20) 3 89 (17) 72 (13) 102 (12) 396 (10) 102 (12) 252 (8) 627 (21) 4 229 (17) 184 (13) 241 (12) 858 (10) 215 (12) 507 (8) 1263 (21) 5 594 (17) 473 (13) 576 (12) 1765 (10) 572 (12) 915 (8) 2399 (22) 6 1426 (18) 1061 (13) 1221 (12) 3345 (10) 1370 (15) 1737 (8) 4044 (22) 37 Table 3.11: Scaled solution time (# of iterations) vs order of polynomial chaos for log-normal random field, advection-diffusion equation, Stoch. dim=4, = 0:1 PC Preconditioners for GMRES GS, Jacobi Solvers order MB AGS AJ GS KP GS Jacobi 2 25 (15) 23 (13) 36 (12) 129 (8) 36 (12) 94 (7) 202 (16) 3 77 (16) 65 (13) 97 (12) 353 (9) 96 (12) 252 (8) 564 (19) 4 229 (17) 184 (13) 241 (12) 858 (10) 215 (12) 507 (8) 1263 (21) 5 733 (18) 585 (14) 645 (12) 1882 (10) 640 (12) 1030 (9) 2631 (24) 6 2157 (19) 1664 (14) 1784 (13) 4085 (10) 1938 (15) 1731 (9) 4821 (26) Table 3.12: Scaled solution time (# of iterations) vs standard deviation () for log-normal random field, advection-diffusion equation, Stoch dim = 4 and PC order = 4 Preconditioners for GMRES GS, Jacobi Solvers MB AGS AJ GS KP GS Jacobi 0.10 229 (17) 184 (13) 241(12) 858 (10) 215 (15) 507 (8) 1263 (21) 0.15 275 (20) 214 (15) 267 (13) 1020 (12) 266 (13) 634 (10) 2083 (34) 0.20 333 (24) 255 (18) 285 (14) 1177 (14) 285 (14) 825 (13) 3889 (62) 0.25 391 (28) 285 (20) 305 (15) 1415 (17) 304 (15) 1018 (16) 17511 (254) 0.30 464 (33) 328 (23) 361 (18) 1656 (20) 362 (18) 1210 (19) Div 0.35 529 (38) 387 (27) 660 (32) 1979 (24) 664 (32) 1464 (23) Div 38 Chapter 4 Basis adaptation - dimension reduction 4.1 Introduction The curse of dimensionality is a significant challenge that remains at the forefront of scientific computing. We address it in the present paper by recognizing the low dimensional structure of many quantities of interest (QoI) even in problems that have been embedded, via parameterization, in very high-dimensional settings. We take advantage of the geometric structure associated with the parameterization to discover very low-dimensional structure around which the probabilistic content of the QoI is concentrated. Recent advances in sensing technologies have enabled querying the natural world with reso- lutions that continue to shed new light on many physical processes that underlie phase transfor- mation and other instabilities. The mathematical formulation of these problems typically involves systems of coupled equations parameterized with stochastic processes that attempt to capture the effect of unmodeled features on the state variables. This stochastic parameterization naturally embeds the problem in a high dimensional parameter space where methods of white-noise cal- culus have recently been adapted to characterizing and estimating the solution of these equations 39 Ghanem and Spanos (1991), Holden et al. (1996), Xiu and Karniadakis (2002). While this setting served to develop general insight into the significance of the stochastic component on the pre- dictability of associated physical processes, it continues to present challenges when dealing with realistic problems that involve very high-dimensional parameterizations. The difficulty is asso- ciated with the observation that, in addition to the fine spatio-temporal discretization required to resolve underlying physical phenomena, resolving the dependence in parameter space increases the number of degrees of freedom exponentially with the number of parameters involved. Some attempts have been made in recent years to address this problem by pursuing sparse Doostan and Iaccarino (2009), Doostan and Owhadi (2011), Schwab and Gittelson (2011) or multilevel Doostan et al. (2007) representations that have succeeded in reducing the computational effort by an order of magnitude. Given the exponential growth with dimension, however, the residual effort remains often prohibitive. This complexity is inherited from the mathematical construction used to describe the physi- cal processes involved, and often belies the simplicity of many quantities of interest (QoI) upon which decisions are eventually based. These QoI’s have typically slow dynamics and pertain to a small number of functionals of the solution. Thus while describing the full solution field as a high-dimensional mathematical object may seem conceptually reasonable, a similar character- ization of inherently low-dimensional quantities seems like an overkill and begs for a rational reduction. Equation-free and projective integration approaches C. Gear et al. (2003) present a useful methodology for tackling these problems when the functional dependence of the QoI on the parameterization of the original problem is not significant to the QoI. In this paper we proceed along different lines, capitalizing on the mathematical structure of Gaussian Hilbert spaces and the associated Homogeneous Chaos space obtained from them through a tensor product construction, 40 to obtain a reduced model that, while capturing the probabilistic content of the QoI, characterizes its functional dependence on the original parameterization. This can be important for carrying out deterministic and stochastic sensitivity analyses and worth of information studies. The approach hinges on Gaussian parameterization of the problem. This does not signify that the parameters in the problem are assumed to be Gaussian, but rather that they are functions of some underlying Gaussian variables or processes. In a first step, an isometry is applied to an underlying Gaussian variables to induce a particular desired structure in the representation of the solution. In a second step, a reduction through projection of the resulting representation is carried out that has a prob- abilistic interpretation, thus providing further justification for the algebraic manipulations. The methodology is demonstrated on an algebraic problem and a problem associated with an elliptic equation from elasticity with stochastic coefficients. 4.2 Polynomial Chaos Decomposition Consider an experiment described by the probability triple ( ;F;P ) and = ( 1 ; ; d ) a set of independent Gausssian random variables on . LetH denote the Gaussian space spanned by andF(H) the-algebra generated byH. For notational simplicity, letL 2 ( ) be a shorthand for L 2 ( ;F(H);P ). Denote byH :n: the space of d-dimensional polynomials in of exact degree n. When generated by Gaussian variables,H :n: is also referred to as the Homogeneous Chaos of ordern. Then it can be shown thatL 2 ( ) = L n H :n: Janson (1999). Given this construction, any basis inH induces a corresponding basis inH :n: andL 2 ( ). The main contribution of this paper is to adapt the hilbertian basis in L 2 ( ) by a judicious basis transformation in the underlying Gaussian space H. Introducing the multi-index = ( 1 ; d ), then each H :n: is spanned 41 byfh ; jj = ng, the set of multivariate Hermite polynomials of order n. We introduce a normalization of these polynomials () =H ()= p !; kH k 2 L 2 ( ) =!; (4.1) and considerq : R d 7! R such that R R d (q(x)) 2 e 1 2 x T x dx <1. Thenq2 L 2 ( ) and can be represented in a Polynomial Chaos decomposition of the form Ghanem and Spanos (1991), q p () = X jjp q (); (4.2) where lim p!1 q p () = q() in the L 2 ( ) norm. In the above,jj = P d i=1 i and ! = Q d i=1 i !. In order to guide our basis adaptation process, we next introduce a few sets of multi-indices. Denote byI p the set of alld-dimensional multi-indices of degree less than or equal top, byE i the subset ofI p corresponding to non-zero indices only at thei th location, and letE = S d i=1 E i . Also letE ij (i6=j) denote the subset ofI p with zeros except at locationsi andj. Moreover, lete i and e ij denote the unit vectors inE i andE ij , respectively, where the non-zero entries are equal to 1. 4.2.1 Change of Basis LetA be an isometry onR d and let be defined as =A AA T =I : (4.3) 42 Also, letq A be the image ofq under this mapping,q A () = q(). Since is also a basis inH, then th order Hermite polynomials in spanH :n . Letting A () = (); (4.4) and given the respective expansions ofq andq A in the forms, q() = X 2Ip q (); q A () = X 2Ip q A (); (4.5) results in q A = X 2Ip q ( ; A ); (4.6) where ( ; ) is the inner product defined as expectation of the product relative to the Gaussian measure. 4.2.2 Reduction via Projection Consider a subspaceV I ofL 2 ( ) spanned byf A ;2Ig, whereII p . The projection ofq A onV I , denoted byq A;I , can be expressed in the form q A;I () = X 2I X 2Ip q ( ; A ) (): (4.7) Alternatively, the projection ofq A onV I can be expressed with respect tof ()g, resulting in q A;I () = X 2Ip q I (): (4.8) 43 This yields, q I = X 2I X 2Ip q ( ; A )( A ; ) (4.9) Denotingq I =fq I ; 2Ig andq =fq ;2I p g, the projectionq I is clearly the restriction ofq toV I . The accuracy ofq I as an approximant toq is thus dependent both on the choice ofA and onI. Clearly, asA approaches the identity map onR d andI approachesI p ,q I approaches q. Also, ifA is fixed as the identity map, then this approximation merely consists of selecting particular terms in the standard polynomial expansion of q. The introduction ofA provides a parameterization of the projections and hence additional freedom for adapting the approximation. Introducing the GramianC defined by its componentsC = ( ; A ), the error associated with a projection onV I can be expressed as, qq I = ICC T q: (4.10) It is clear that the null space ofICC T and the spectrum ofC are shaped by bothA andI. Our construction does not attempt to minimize any measure of this error, but rather seeks probabilistic arguments for the construction ofA and the selection of the projection setI. 4.3 Special Cases Two special cases are now considered whereA is used to introduce knowledge aboutq into the process of basis adaptation. 44 4.3.1 Case where Gaussian components of q are known If the projection ofq on the Gaussian subspace,H, is known, possibly through some extraneous numerical effort, a very interesting reduction ensues. Specifically, letA be constructed such that 1 = X 2I 1 q () = d X i=1 q e i i (4.11) so that the first component of captures the complete Gaussian component ofq. Let the remain- ing components of the vector be constructed through a Gram-Schmidt procedure, yielding the matrixA. Furthermore, let the setI introduced previously consist of all multi-indices that refer- ence only 1 up to degreep, and is thus specified asI =E 1 . In this case,q A () andq A;I () are given respectively as, q A () =q A 0 +q A e 1 1 + X 1<jjp q A () =q A 0 +q A e 1 1 + X 1<jjp 2E 1 q A () + X 1<jjp 62E 1 q A (): (4.12) and q A;I () =q A 0 + X 2E 1 q A () (4.13) and the error in the representation is given as, q A ()q A;I () = X 1<jjp 62E 1 q A () (4.14) 45 which has mean zero and is orthogonal to all the Gaussian terms. Thus, starting with a knowledge of the Gaussian component of q (i.e. its projection on the Gaussian chaos), an isomeryA is constructed such that this whole projection is carried by a single dimension. The choice ofI is then merely adapted to the situation where the probability measure ofq is concentrated around this leading dimension. Clearly, the setI can be enlarged by activating successive’s in the list of multi-indices. In the limit when all are active,I =I p and an exact representation is obtained. 4.3.2 Case where quadratic components of q are known It is assumed in this subsection that the quadratic components ofq are known. An expression of q that highlights the second order terms takes the form, q() =q 0 + d X i=1 ^ q i i + d X i=1 d X j=1 ^ q ij ( i j ij ) + X 3p q () (4.15) where ^ q i = q e i , ^ q ij = q e ij = p 2, and ^ q ii = q 2e i = p 2. The second order truncation of equation (4.15) can be re-written as q II () = (q 0 d X i=1 ^ q ii ) + d X i=1 ^ q i i + d X i=1 d X j=1 ^ q ij i j : (4.16) Let matrixS have entries ^ q ij , and let =A where ASA =D (4.17) 46 whereD is a diagonal matrix whose elements are denoted asd i . ThusA andD are, respectively, the matrix of eigenvectors and the eigenvalues of the symmetric matrixS, and thusA is clearly a unitary operator. Noting that the trace ofA is equal to the trace ofD,q A () takes the form, q A () =q 0 + d X i=1 b i i + d X i=1 d i ( 2 i 1) + X jj>2 q (); (4.18) whereb i = P d j=1 a ij ^ q j anda ij denotes elements ofA. LettingI =E, the projection onV I takes the form, q A;I =q 0 + d X i=1 X 2E i q A (); (4.19) which captures exactly the second order effects while including the higher order effects only along one-dimensional subspaces. This approximation can be achieved usingd one-dimensional poly- nomial approximations, each of orderp. Clearly, constructing the operatorA requires additional work to evaluate the second order approximation of function q. The error associated with this approximation can be expressed as q A ()q A;I () = X 2<jjp 62E q A (): (4.20) The error reflects the interaction between the i components with polynomial order greater than 2. It can be reduced by augmentingV I with these higher-order interactions. 47 4.4 Application to Equations with Random Coefficients Consider the situation whereu2V is the solution of an equation with random coefficients,V is some appropriate function space that has at least the structure of a topological vector space. Let L denote the operator representing this equation and which is parameterized by2R d , a vector of independent Gaussian variables. We note that this does not mean that the parameters of the underlying model are themselves Gaussian variables, but merely that they have been described as functions of such variables. We are interested in characterizing the probabilistic structure, through its chaos decomposition, of a functional,q, onV . The problem could then be posed as characterizing the solution of the following equation subject to sufficient additional constraints (such as boundary or initial conditions), L [u] = 0 8u2V L 2 ( ) q() = f(u()); (4.21) where the dependence ofu andq on is emphasized. We assume throughout thatL andf are such thatu andq are square-integrable with respect to the density of, and they therefore admit polynomial chaos decompositions. 4.4.1 Non-Intrusive Implementation A good transformationA is such that the functionq does not vary significantly except along a small subset of in which case integrals ofq over the whole space can be expressed as much lower-dimensional integrals with respect to that subset. 48 The coefficients in the chaos decomposition ofq can be expressed as an integral which is then approximated through quadrature as, q = Z R d q() ()e 1 2 T d s X r=1 q( (r) )w r ; 2I p (4.22) wheref (r) g andfw r g are the coordinates and associated weights for ad-dimensional quadrature rule withs points, and the integrand is the solution of the following deterministic problem L (r)[u] = 0 u2V; r = 1; ;s q( (r) ) = f(u( (r) )); r = 1; ;s: (4.23) The computational burden of carrying out the above high-dimensional quadrature increases with jj and quickly becomes prohibitive. This procedure can be used, however, to construct a first or- der (jj = 1) or second order (jj = 2) representation, from which an isometryA is constructed as described above. Assuming, for instance, that the projection of the solution on the Gaussian subspace has been computed using ad-dimensional quadrature as described in this section, equa- tion (4.11) is used to constructA, and then a one dimensional representation ofq is sought in the form, q =q 0 + X 2E 1 q (): (4.24) The coefficients q in this expansion can now be obtained through one-dimensional s o -point quadrature with respect to 1 in the form, 49 q = so X r=1 q( (r) ) ( (r) )w r ; 2E 1 ; (4.25) where (r) = ( (r) 1 ; 0; ; 0), andq( (r) ) is obtained by solving the following equations, (r) = A 1 (r) ; r = 1; ;s o L (r)[u] = 0 u2V; r = 1; ;s o q( (r) ) = f(u( (r) )); r = 1; ;s o : (4.26) In this case, and as previously explained,A is merely constrained by requiring = A, and its construction completed through an arbitrary Gram-Schmidt procedure. Clearly, the non- uniqueness ofA will be manifested in the first of equations (4.26). For quantities of interest that are strongly dominated by 1 , however, it is expected that this arbitrariness in specifying is inconsequential in evaluating the integral. 4.4.2 Intrusive Implementation In some instances, reduction of the multi-dimensional integral as indicated above may not be justified. In that case, a procedure as explained in this subsection can be implemented. It will again be assumed that some effort has already been expanded in computing the isometryA that is adapted toq, together with an approximating subspaceV I . Without loss of generality, and to keep the focus of the presentation on the basis adaptation procedure, it is assumed throughout this section thatL is a linear operator. Describingu in its polynomial chaos decomposition, and 50 requiring the error to be orthogonal to the approximating space results in the following system of deterministic equations to be solved for the coefficientsu , X 2Ip (L [u ] ; ) = 0 8 2I p : (4.27) where (:;:) is the inner product introduced in equation (4.6). In some instances, the operatorL lends itself to a chaos decomposition, whereby we can introduce the deterministic operatorsL in the form, L = (L ; ); 2I p ; (4.28) resulting in X 2Ip X 2Ip L [u ] ( ; ) = 0 8 2I p : (4.29) If the solution is instead represented relative to the adapted basis the following set of equa- tions ensues, X 2Ip X 2Ip L [u ] A ; A = 0 8 2I p ; (4.30) where A ; A is an integral with respect to the Gaussian measure, and represents an intrin- sic property of Hermite polynomials that could be evaluated ahead of the numerical solution of the equations. While analytical expressions for this integral can be obtained, it may be most effective to tabulate it, for different values ofA, using numerical quadrature inR d . The integrand in this quadrature is an inexpensive function whose value is independent of the physics of the problem or the accuracy of the discretization scheme of the space-time continuum (onceA is known). 51 4.4.3 Adapted Representation of Model Parameters It is reasonable to expect that different quantities of interest should poll the parameters of the model differently for relevant information. Given the mathematical formalism adopted in this paper for model reduction, namely basis adaptation and orthogonal projection, this issue can be stated clearly and concisely. For that, we revisit the intrusive implementation in section (4.4.4.4.2) and we express the random operatorL in terms of its polynomial chaos decomposition relative to, with equations (4.28) and (4.30) replaced by the following equations, respectively, L = (L ; A ) 2I; (4.31) and X 2I X 2Ip L [u ] A A ; A = 0 8 2I; (4.32) where we have, in addition, accounted for the projection of the solution onV I . Since both and are standard orthogonal gaussian vectors, A A ; A = ( ; ) ; (4.33) and equation (4.32) can be rewritten as, X 2I X 2I L [u ] ( ; ) + X 2I X 62I L [u ] ( ; ) = 0 8 2I: (4.34) 52 If the multi-index setI is chosen to be orthgonal to its complement inI p (which is indeed the case for the two proposed choices of the cases consi dered in sections 4.3.4.3.1 and 4.3.4.3.2, then the second double summation in the last equation vanishes, and one is left with a problem identical in structure to the original problem stated in equation (4.29). Thus, if we had the foresight of knowing the subsetI before solving the problem, and expanding the model parameters along that multi-index set, the solution of that problem for the QoI would coincide with the solution for the QoI obtained from a full-dimensional solution of the problem. In general, however, it is not possible to expand the model parameters, ahead of time, with respect to the multi-index setI since this set characterizes the yet to be determined solution. The value of the foregoing analysis, however, is to recognize that to any given QoI adapted to the multi-index setI, there exists a commensurate characterization of the model parameters, obtained via projection, that quantifies how much of these properties is needed for the accurate estimation of the QoI. 4.5 Numerical Example We now consider two examples, demonstrating the applicability of the foregoing analysis to prob- lems specified, respectively, by algebraic and partial differential equations. 4.5.1 Algebraic Equation Consider the function ^ h :R d 7!R given by the expression, ^ h( ^ 1 ; ; ^ d ) = d 1 +b d Y i=1 (1 +a i ^ i +b i ^ 2 i ) d X j=1 d Y i=1 i6=j (1 +a i ^ i +b i ^ 2 i ) (4.35) 53 Figure 4.1: Probability density function ofh(). Results using a full 10-dimensional solution and the reduced approximations using prior knowledge on the Gaussian content and quadratic content of QoI. which represents the effective modulus ford springs in series, each of modulus (1 +a i ^ i +b i ^ 2 i ). Correlation is introduced between theN (0; 1) Gaussian random variables ^ i in the form of Ef ^ i ^ j g = 2 e jijj=l . A Karhunen-Loeve expansion is first affected to transform the setf ^ i g into a set of independent random variablesf i g. Alld terms are retained in the resulting expan- sion, and the induced mapping,h, has anL 2 expansion with respect to the Hermite polynomials with coefficientsfh g. Evaluating the coefficients in the full third order 10-dimensional expan- sion required a level 5 Smolyak quadrature (8761 function evaluations). First and second order expansions are obtained by specifying, respectively,I p =I 1 andI p =I 2 . Coefficients in these expansions can be obtained via level 2 (21 function evaluations) and level 3 (221 function eval- uations) Smolyak quadratures, respectively. MatrixA is then obtained as indicated by equations (4.11) and (4.17), and the setI is identified, respectively, withE 1 andE. The numerical results below are associated with statistically identical springs withb i =b = 1,a i = 0:5, = 0:05 and 54 p = 3. The correlation length of the springs was taken asl = 4 We note that the effective modulus becomes increasingly Gaussian as the correlation length approaches zero. Figure (4.1) shows the results for the probability density function ofh, for both linear and quadratic prior knowledge. The quality of the approximation in the tail region is also shown. 4.5.2 Linear Elasticity with Nonlinear Quantities of Interest As a second example, consider an elastic material occupying a square domain and subjected to a uniform pressure along its boundary with an empty circular inclusion. Due to symmetry, only a quarter of the domain is analysed with suitable boundary conditions, as depicted in Figure (2). The material behavior is assumed to be characterized in the realm of linear elasticity as a plane stress problem. It is further assumed that the modulus of elasticity of the medium is a spatially-varying stochastic process while Poisson ratio is a deterministic constant over space. This assumption is clearly inconsistent with linear elasticity since the modulus of elasticity and Poisson ratio are both dependent on the same microstructure. Our working assumption, however, is consistent with observations of much smaller statistical scatter in Poisson ratio than the modulus of elasticity, and can be viewed as resulting in a model that is a perturbation to mathematical elasticity. The governing equations are discretized using the finite element method according to the mesh shown in Figure (4.2). The modulus of elasticity is assumed to be the exponential of a Gaussian field and is defined by its normalized spatially uniform mean value,E 0 = 1N=m 2 :, an isotropic exponentially decaying covariance function with correlation coefficient equal toL=8 where L = 1m is the length of one side of the square domain, and a coefficient of variation equal to 0:5. Poisson ratio is taken to be 0:25, and the pressure along the boundary with the inclusion is taken as 0:1N=m 2 . A Karhunen-Loeve expansion is used to discretize the Gaussian 55 Figure 4.2: Elastic domain with random Young’s modulus and subjected to uniform internal pres- sure. PointA (0.2977,0.1417) is highlighted. process associated with the modulus of elasticity, with 10 terms retained in its expansion (d = 10). The exponential of this representation is then developed in a 3 rd order 10-dimensional Hermite polynomial expansion. The solution to the elasticity equations with this representation of material properties consists of two displacement fields (one along each spatial dimension) that form a 2- variate spatial stochastic process. A polynomial chaos expansion of this process up to order 3 in all 10 dimensions is carried out using level-5 Smolyak quadrature, as a benchmark against which to verify the model reduction results. This provides our estimate for the processu. We then focus attention on two separate scalar quantities of interest (QoI). The first QoI is a linear functional ofu and consists of the displacement in the horizontal direction at one point in the domain. The second QoI is a nonlinear functional and consists of the maximum vonMises stress over the two dimensional spatial extent of the problem. V onMises stress,, at any given point in the material is given as = p 2 1 1 2 + 2 2 where 1 and 2 are the two principal stresses at that location. As a motivation to the reduction methodology developed in this paper, we computed the Karhunen-Loeve expansion of the Gaussian component of thex-displacement process, and used the Gaussian random variable associated with its dominant mode to construct a 1-dimensional 56 polynomial chaos approximation of the first QoI described above. Figure (4.3) compares the results from this 1-dimensional approximation against the full 10-dimensional evaluation of the QoI and the approximations according to the two methods described in this paper. It is clear that the reduction using the dominant KL variable fails to capture the probabilistic measure of the QoI, emphasizing the need for more rational adaptation strategies. The probability density functions of the first QoI obtained using the two procedures described earlier in the paper are shown in the same figure. The linear (Gaussian) approximation is obtained using a level 2 Smolyak quadrature whereas the quadratic approximation required a level 3 rule. The match is amazing even in the tail of the distributions, using just the one-dimensional approximation. Figure (4.4) shows a comparison of the polynomial chaos coefficients associated with a 10-dimensional approximation and the approximation obtained using a sequence of 10 one-dimensional solutions adapted to the quadratic components of QoI. It is clear that the terms being neglected in the proposed 1d reduction have minor contribution to theL 2 error. Figure (4.5) shows analogous results for the second, nonlinear, QoI. In this case, an approx- imation built around the Gaussian content of the QoI is insufficient to completely characterize its probability density function. The reduced model obtained by including information from the second order terms, however, provides an excellent approximation. The analysis presented in section 4.4.3 is applied to the present problem. Figure (4.6) shows the the first six polynomial chaos coefficients of Young’s modulus associated with an expansion of the lognormal process. Note the global features in these coefficients. Figure (4.7) shows the same quantities for the problem adapted to the QoI specified by the vonMises stresses at pointA. Note the localization around pointA in the polynomial chaos coefficients, induced by the quantity of interest. 57 Figure 4.3: Probability density function ofx–displacement at pointA. Results using a 1) full 10- dimensional solution and the approximations using the 2) dominant component in the Karhunen- Loeve expansion of the stochastic process associated with the solution, 3) the 1-d approximation associated with Gaussian prior knowledge, and 4) sequence of 10 1-d approximations associated with quadratic prior knowledge. Figure 4.4: Coefficients of the PC expansion of thex-displacement at pointA. Results using a full 10-dimensional solution and the reduced approximation using prior knowledge on the quadratic content of QoI. 58 Figure 4.5: Probability density function of the maximum vonMises stress over the spatial do- main of the problem. Results using 1) full 10-dimensional solution and the approximations using the 2) 1-d approximation associated with Gaussian prior knowledge, and 3) sequence of 10 1-d approximations associated with quadratic prior knowledge. This correspondance between a QoI and certain features of model parameters should have implications on the characterization of material properties from experimental data and could serve as a guide for the development of adapted regularization strategies for the associated inverse problem. That invesigation, however, is beyond the scope of the present work where we limit ourselves at noting the presence of such an interplay. 4.5.3 Random eigenvalue problem Eigenvalue analysis has very important application in structural dynamics to find fundamental frequency of the dynamical system and and corresponding eigenmodes. When parameters of the system of the system are random, the eigenvalue problem becomes, random eigenvalue analysis. 59 Figure 4.6: The first six coefficients of the elasticity modulus with respect to a polynomial chaos expansion in the original gaussian variables. Figure 4.7: The first six coefficients of the elasticity modulus with respect to a polynomial chaos expansion in the one-dimensional gaussian variable 1 . 60 Deterministic eigenvalue problem Let K 2 R nn and M 2 R nn be stiffness and mass matrices obtained from finite element discretization. To find fundamental frequency and corresponding mode shapes of the dynamical system, following eigenvalue problem has to be solved, K = M ; 2R; 2R n ; (4.36) where, and are eigenvalue and eigenvector respectively. Random eigenvalue problem When there is uncertainty in the system, then the eigenvalue problem becomes following random eigenvalue problem, K()() =()M()(); 2 where,K()2R nn ,M()2R nn are random stiffness and mass matrices,()2R;()2 R n are random eigenvalue and corresponding random eigenvector and ( ;F;) is underlying probability space with sample space ,-algebra,F and probability measure such that, each sample point2 . Polynomial chaos approach Random eigenvalue problem discussed in section 4.5.3 can be solved using polynomial chaos expansion method discussed in section 2. Stiffness matrixK, mass matrixM, eigenvalue and eigenvector can be expanded in terms of polynomial chaos inL 2 ( ;F;) asK = P P1 i=0 K i i (), M = P P1 i=0 M i i (), = P P1 i=0 i i () and = P P1 i=0 i i (), where =f i ; ; d g T . 61 Random eigenvalue problem can be solved using stochastic Galerkin or stochastic collocation methods. Stochastic Galerkin projection of random eigenvalue problem results in following cou- pled system of equations, P1 X i=0 P1 X j=0 Ef i j k gK i j = P1 X i=0 P1 X j=0 Ef i j l k g i M j l ; k = 0; ;P 1: (4.37) Stochastic collocation method involves solving deterministic eiegenvalue problems at a given set of quadrature points. K( (q) )( (q) ) =( (q) )M( (q) )( (q) ); q = 1; ;q n In this work stochastic collocation methods is used for testing basis adaptation methods for ran- dom eigenvalue problem. The basis adaptation methodology can be readily applied to random Figure 4.8: 2-d cantilever beam meshed with quad4 elements, is fixed on left end and eigenvalue analysis is performed eigenvalue problem to represent the eigenvalue in one random variable and solve a one dimen- sional random eigenvalue problem to compute the eigenvalue accurately. These approaches are tested in two structural systems, a two dimensional cantilever beam shown in figure 4.8 and a wind turbine blade sgown in figure 4.9. Figure 4.10 shows the probability density functions of first six 62 Figure 4.9: Wind turbine blade meshed with quad4 elements, is fixed in all 6 dof on one end and eigenvalue analysis is performed eigenvalues of a two dimensional cantilever beam and figure 4.12 shows the pdfs of first six eigen- values of the wind turbine. It can be observed from both cantilever beam and the wind turbine blade problems that the eigenvalues can be captured very accurately using one dimensional basis adaptation. Although the eigenvector corresponding to each eigenvalue can be computed after the eigenvalue is computed, gain in the computational effort has to be investigated. To adapt the basis adaptation to eigenvectos, letf i g r i=1 be a set of deterministic eigenvectors, then a random eigenvector (l) () cane be can be projected on to the deterministic set of nominal eigenvectore f i g r i=1 to get coefficients, (l) () = P1 X i=0 (l) i () (l) () = r X j=1 (l) j P1 X i=0 a (l) ij i (); (4.38) 63 (a) first eigenvalue (b) second eigenvalue (c) third eigenvalue (d) fourth eigenvalue (e) fifth eigenvalue (f) sixth eigenvalue Figure 4.10: pdf of eigenvalues of 2-d cantilever beam in 6d and in 1d then, the coefficients, (l) j can be computed as quantities of interest using basis adaptation meth- ods. The first coefficient (l) 1 , polynomial chaos coefficients computed in reduced dimension form basis adaptation is in very good agreement with those computed in high dimension. How- ever, there is a significant error in other coefficients. It could be due to the lower magnitude of the second and third coefficients. 64 (a) first eigenvalue (b) second eigenvalue (c) third eigenvalue (d) fourth eigenvalue (e) fifth eigenvalue (f) sixth eigenvalue Figure 4.11: pdf of eigenvalues of wind turbine blade in 6d and in 1d (a) (l) 1 (b) (l) 2 (c) (l) 3 Figure 4.12: Polynomial chaos coefficients of of (l) j in 6d and in 1d 65 Chapter 5 Stochastic upscaling 5.1 Introduction In this chapter, stochastic upscaling methodology for a steady state laminar flow through an ar- ray circular solid discs and heat transfer explored. The fine scale problem is modeled using Navier-Stokes equations and solid discs act as heat source (sink). The uncertainty in the system is modeled by treating thermal conductivity of the solid discs as a random filed. At coarse scale the medium is modeled as a porous medium and Darcy-Brinkman equations are used. Due to temperature dependent fluid density and fluid viscosity, the uncertainty in the fine scale is prop- agated to fluid flow and parameters such as permeability and thermal conductivity at the coarse scale. V olume average method is adopted to compute the permeability and thermal conductivity at the coarse scale from fine scale solution. Polynomial chaos expansion methodology is used to solve stochastic partial differential equations at fine and coarse scales. Preconditioners devel- oped in chapter 3 are used to solved stochastic Galerkin system of equations efficiently. Basis adaption methods developed in chapter 4 used to further reduce the stochastic dimension while 66 computing the coarse scale parameters, namely, for stochastic dimension reduction for computing permeability random field and thermal conductivity random field. 5.2 Problem statement Many physical problems can be described with coupled stochastic nonlinear partial differential equations. One such problem of interest is fluid flow through porous media coupled with heat equation. At fine scale the porous media problem can be solved using Navier-Stokes equation. However, it is computationally very expensive to solve using Navier-Stokes equations with ran- dom parameters. For this reason, the system upscaled to porous media equations and parameters of the porous media are computed from fine scale solutions. Uncertainty at the fines scale is modeled by treating the thermal conductivity of the solid discs as a log-normal random field. The uncertainty at the fine sale is propagated to the coarse scale which can be modeled by treating the permeability and thermal conductivity as random parameters. To solve these stochastic partial dif- ferential equations at fine and coarse scales, use spectral stochastic finite element method Ghanem and Spanos (1991) is used. Stabilized finite element methods Tezduyar (1992) are used to circum- vent the problem of instability due to spurious pressure modes that arise due to equal order in- terpolation for pressure and velocity fields in finite element discretization of incompressible fluid flow equations. Figure 5.1 shows fine scale and coarse scale domain modeled in 2-dimensions. 67 Figure 5.1: Fine scale domain modeled with an array of solid discs and coarse scale domain modeled as porous media 5.2.1 Navier-Stokes equations at fine-scale At fine scale, two dimensional system is modeled as an incompressible flow though an array of solid circular discs. Heat transfer between solid circular discs and the fluid is modeled through a coupled heat equation. As the fluid flows through the space between the discs, heat is trans- ferred from discs to the fluid. Uncertainty due to thermal conductivity of the solid discs is propagated to the fluid velocity, pressure and temperature through, temperature dependent fluid viscosity and fluid density. This physics at the fine scale can be modeled with Navier-Stokes equations Nakayama et al. (2002) and heat equation as f (ru) = 0 (5.1) f (ru)u =rP + f r 2 u + f f (5.2) f c p;f [r (uT )] =k f r 2 T (5.3) 68 where, equations 5.1, 5.2 and 5.3 are the continuity, conservation of momentum (Navier-Stokes) and energy equations respectively. In these equations,u is the velocity vector of the fluid,P is pressure in the fluid,T is the temperature andf is the forcing term. 5.2.2 Darcy-Brinkman equations at coarse-scale At coarse scale the domain is treated as a porous medium with heat source in the solid region and the physics is modeled with the Darcy-Forchheimer equation Schmidt et al. (2010). The equations governing flow through porous media are as follows, f (r ~ u) = 0 (5.4) f f (r ~ u)~ u =r ~ P + f r 2 ~ u f K ~ u + f ~ f (5.5) r C~ u T = k eff r 2 T + _ Q (5.6) where, equations 5.4 and 5.5 are the continuity and momentum equations and 5.6 is the heat equation under thermal equilibrium between fluid and solid phases. where, ~ u is the velocity of the fluid, ~ P is the pressure in the fluid, is the porosity of the medium, C is the specific heatK is the permeability, k eff is effective heat conductivity of porous media and _ Q s is heat source in the solid. 69 At the coarse scale, uncertainty in the system can be represented by treating the permeability and thermal conductivity as random parameters. To solve these stochastic partial differential equations at fine and coarse scales, spectral stochastic Galerkin projection method Ghanem and Spanos (1991) will be used. Stabilized finite element methods Shadid et al. (2006), Tezduyar (1992) are used to circumvent the problem of spurious pressure modes that arises due to equal order interpolation in solving incompressible Navier-Stokes equations. 5.3 Stochastic upscaling Stochastic upscaling is achieved by computing the coarse scale random parameters using volume average method from the fine scale solution. Further, stochastic basis adaptation methods are used to represent the coarse scale parameters in reduced stochastic dimension. 5.3.1 Permeability The coarse scale permeability is computed by volume average method Nakayama et al. (2002). At low velocity the Darcy law can be used to compute the permeability @hPi @x i =k 1 ij hu j i (5.7) Where,hPi andhui are volume averaged pressure and velocity obtained from fine scale solution. V olume average is computed as hai = 1 V Z V adV (5.8) 70 where V is the representative volume considered for averaging. The solution of the fine scale stochastic partial differential equations are obtained using non-intrusive Smolyak quadrature method. At each quadrature point, a deterministic fine scale problems is solved and fine scale, velocity, pressure and temperature distribution are obtained. From the velocity and pressure solution, a permeability filed is constructed for each solution. Due to very small variability of pressure gra- dient is y-direction, the permeability is considered only along x-direction. From the quadrature points, polynomial chaos expansion of the permeability filed and higher order statistics can be readily computed. Figure 5.8, shows the mean and standard deviation of the permeability field. Thus the uncertainty due to thermal conductivity in the fine scale problem is propagated to the coarse scale problem through the coarse scale permeability parameter. Further, the covariance function of the permeability field is constructed and Karhunen-Lo` eve analysis is performed to find the dominant eigenvalues and eigenmodes. 5.3.2 Thermal conductivity Thermal conductivity of the porous media is upscaled from the finescale heat equations. At coarse sclae thermal conductivity relates the Darcy velocity to the temperature distribution at the coarse scale. The coarse scale heat equation is written as r C~ u T = k eff r 2 T (5.9) Above relation can be rewritten as Cr ~ u T = k eff r (r T ) (5.10) 71 Where, and C density of the fluid and specific heat at the coarse scale. k eff is effective thermal condictivity of the porous media ~ u is the coarse scale velocity and T is the coarse scale tempera- ture. Writing above equation interms of fine scale solution,u andT , C Z D r (uT )dA = k eff Z D r (rT )dA: (5.11) Using Gauss-divergence theorem, area integrals in the above equations can be converted intoline integrals as follows, C Z C ~ n (uT )ds = k eff Z C ~ n (rT )ds; (5.12) where,u andT are velocity and temperature from the fine scale solution. 5.4 Basis adaptation to find QoI In this section, we discus the basis adaptation methods to compute the permeability random filed in reduced dimension. From the volume average method discussed in section 5.3, permeability random fieldk(x;) can be computed from the fine scale solution as following polynomial chaos expansion, k(x;) = p X i=0 k i (x) (); (5.13) 72 where, =f 1 ; ; d g T is the input random vector obtained from truncated Karhunen-Loeve expansion of the thermal conductivity random field at the fine-scale. Thermal conductivity ran- dom field is modeled with log-normal distribution and hence, it can be expanded in terms of Hermite polynomials of standard Gaussian random variables. 5.4.1 Linear basis adaptation Linear basis adaptation for permeability random field is pursued by computing a new set of ran- dom variables = A using an isometryA such thatAA T = I. A is not unique and one way to constructA is using the Gaussian components in the PC expansion of the quantity of interest (QoI). In the current problem, the permeabilityk(x; at a spatial locationx j can be computed as a QoI in the reduced dimension. That is , 1 = d X i=1 k i (x j ) i ; (5.14) here,k i (x j ), after normalization forms the first row of matrixA. Rest of the rows in the matrix A can be constructed using Gram-Schmidt process. Now the permeability at location,x j can be constructed in reduced stochastic dimension, with very few randm variables. However, using linear basis adaptation, new set of has to be computed for each spatial location. 5.4.2 Quadratic basis adaptation Quadratic basis adaptation procedure can be pursued when, the quadratic terms in the polynomial chaos expansion QoI are known. In this case the isometry matrixA is obtained as follows, ASA T =D; (5.15) 73 , where D is the diagonal matrix, whose elements d i are eigenvalues and A is the matrix of eigen vectors of symmetric matrix S. Here, matrix S contains the quadratic terms in the PC expansion of the QoI. This transformation eliminates, the cross terms in the PC expansion and hence reduces the d-dimensional stochastic Galerkin system in to d one dimensional systems in i ;i = 1; ;d. As seen in the case, linear basis adaptation, quadratic adaptation also, uses new set of at each location and hence, it is difficult to construct entire permeability random field simultaneously. To address this issue, we pursue, basis adaptation using, KL expansion of the Gaussian part of the random filed which is discussed in the next section. Figures 5.10 show the permeability at points 3 and 20, computed using linear and quadratic basis adaptation separately. 5.4.3 Basis adaptation with KL expansion When, the QoI is a scalar, it can be effectively represented in one random variable, using linear or quadratic basis adaptations discussed in the previous sections. However, if the QoI is a random field, above methods can’t be used to simultaneously construct entire random filed. For this reason we use, when the Gaussian part of the permeability random fieldk g = k 0 (x) + P d i=1 k i (x) is known, we construct the isometry from KL expansion of the covariance cov kg of the Gaussian part of the permeability field, cov kg (x 1 ;x 2 ) = d X i=1 k i (x 1 )k i (x 2 ); (5.16) and k g (x;) = ~ k 0 + d X i=1 p i ~ k i (x) i ; (5.17) 74 where, i are the new set of standard normal random variables and can be computed as i = 1 p i Z k g (x;) ~ k i (x)dx; (5.18) Figure 5.9(b) shows the decay, of the eigenvalues ( i ). It can be observed that the decay of random variable is very sharp and very few number (2-3) of random variables is sufficient to represent the entire permeability random field. Figure 5.12 shows the pdf of permeability computed at a point using, full ten dimensional, reduced one, two and three dimensional. It can be observed that, three random, variables 1 , 2 and 3 are sufficient to capture the pdf accurately. 5.5 Results and discussion A two dimensional fluid flow problem is considered with the geometry shown in figure 5.2. Red circles are the solid discs with heat source and the rest of the domain is fluid flow domain. Ther- mal conductivity of the solid discs is modeled as a log normal random field with mean 150 and coefficient of variation 1.0. Temperature dependent fluid density and viscosity are used to main- tain a coupling between fluid flow and heat equations. A Neumann boundary condition of ve- locity, 50cm=s is applied on the left side boundary and velocities at top and bottom boundaries are provided as zero Dirichlet boundary condition. Heat source in the solid discs is provided by specifying a Dirichlet boundary condition at a small region inside the solid disc. Figure 5.3 shows two realizations of the thermal conductivity. Figure 5.9(a) shows the decay of eigenvalues from Karhunen-Lo` eve analysis of thermal conductivity random filed and figure 5.9(b) shows the decay of eigenvalues from Karhunen-Lo` eve analysis of the Gaussian part of the permeability random filed at the coarse scale. Figures 5.4, 5.6 and 5.7 show the mean and standard deviation of the 75 velocity in x-direction, pressure and temperature respectively. Figures 5.8(a) and 5.8(b) show the mean and standard deviation of the coarse scale permeability. Figure 5.2: Fluid and solid regions meshed with 2-d quad4 elements (a) Realization 1 (b) Realization 2 Figure 5.3: Realizations of the thermal conductivity in the circular discs Figure 5.10 shows the pdfs of permeability of the coarse scale parameter at two spatial lo- cations (at point number 3 and 20) computed using 10-dimensional -space and reduced 1- dimensional and 2-dimensional -space. It can be observed the the coarse scale permeability 76 (a) Mean ofux (b) std. dev. ofux Figure 5.4: Mean and std. dev of velocity in x-direction in 10d and order 3 (a) Mean ofuy (b) std. dev. ofuy Figure 5.5: Mean and std. dev of velocity in y-direction in 10d and order 3 (a) Mean of pressure (b) std. dev. of presure Figure 5.6: Mean and std. dev of pressure in 10d and order 3 (a) Mean of temperature (b) std. dev. of temperature Figure 5.7: Mean and std. dev of temperature in 10d and order 3 parameter represented in reduced dimension is very accurate. However, with this method, a new set of has to be identified for each spatial locations. Alternately, figure 5.12 shows the pdfs of permeability computed in reduced dimension using basis adaptation based on Karhunen-Lo` eve 77 (a) mean ofk11 (b) std. dev. ofk11 Figure 5.8: Mean and std. dev. of permeability componentk 11 using volume averaging (a) eigenvalues of fine-scale random field (b) eigenvalues of coarse-scale random field Figure 5.9: Eigenevalues obianed from KL expansion of thermeal conductivity at finescale and the eigen values obtained from KL analysis of Gaussian part of permeability at coarse scale expansion as discussed in section 5.4.3. Using this method, entire permeability is computed simultaneously, but needs 3-dimensional space with 1 ; 2 and 3 . Figure 5.12 shows thermal conductivity computed at point 20 using the basis adaptation method based on Karhunen-Lo` eve expansion. 78 (a) Permeability at point 3 (b) Permeability at point 20 Figure 5.10: Permeability computed at points 3 and 20, using linear and quadratic basis adaptation methods. Here, the basis adapted for each spatial point (QoI) separately (a) Permeability at point 3 (b) Permeability at point 20 Figure 5.11: Permeability computed at points 3 and 20 using basis adaptation method, where new random variables are obtained from Karhunen-Loeve analysis of the Gaussian part of the random filed. Here, basis is adapted for all spatial points (QoI) simultaneously. 79 Figure 5.12: Thermal conductivity computed at points 20 using basis adaptation method, where new random variables are obtained from Karhunen-Loeve analysis of the Gaussian part of the random filed. Here, basis is adapted for all spatial points (QoI) simultaneously. 80 Chapter 6 High resolution micrograph synthesis 6.1 Introduction Recent capabilities in the multi-scale modeling and simulation of material behavior have been matched by technologies for sensing and characterizing materials at distinct spatial scales. Several technical challenges ensue in the process of pairing multi-scale models to multi-scale data, one of which is addressed in this paper namely that of characterizing knowledge about multi-scale material structure in a manner that is conducive to probabilistic conditioning and updating. In this work experimental data provided through micrographs are considered, which are dig- ital images of the microstructure taken through a microscope, showing a magnified view of the material domain, and allowing digital representation of fine-scale features. In this paper, the terms “image”, “texture” and “micrograph” are used interchangeably. High resolution microstructural information is required as a prerequisite for the multi-scale modeling and analysis of compos- ite structures, and in particular to track the nucleation of instabilities in material behavior. As it remains expensive to experimentally obtain high resolution micrographs throughout a computa- tional domain where coarse-scale behavior is being assessed, procedures for gleaning equivalent 81 information from lower resolution images have demonstrated their value in computational mate- rials science Ghosh et al. (2006), Valiveti and Ghosh (2007). This paper contributes to that body of knowledge by developing a statistical procedure that updates low resolution priors with spa- tially localized high resolution data using a parametric texture model and a particle filter. Three approaches are commonly used, in light of their simplicity, for generating high resolution images from lower resolution ones. These are based on nearest neighbor, bilinear and cubic interpola- tion methods Lehmann et al. (1999). In the nearest neighbor interpolation method, a piecewise constant interpolation maps the low resolution to the high resolution image, typically producing blocky images. Bilinear and cubic interpolation methods smoothen these discontinuities, while still allowing for a blur in the high resolution image. Methods based on wavelet based interpo- lation with gradient based enhancement (WIGE) and higher-order polynomial interpolation with gradient based enhancement (PIGE) techniques have been proposed and successfully adapted to problems of computational multi-scale material modeling Ghosh et al. (2006), Valiveti and Ghosh (2007). In WIGE, high resolution images are obtained by interpolating low resolution images using wavelet basis with the interpolated images enhanced using gradient based meth- ods calibrated from a few high resolution images. In PIGE, higher order polynomial bases are used for interpolation. In these methods, gradient based enhancement accounts for the position of calibration images relative to the spatial location of the simulated image. Both WIGE and PIGE methods are shown to be very effective in synthesizing high resolution images from low resolution ones with the aid of few calibrating high resolution images. In Singh et al. (2006), digital micrographs are simulated by representing the microstructures with two-point correlation function. Characterization and simulation of three-dimensional microstructures is also discussed in Groeber et al. (2008a,b), Rollett et al. (2007), Tewari et al. (2004). In these approaches, the 82 synthesized micrographs are deterministic and do not account for uncertainty in experimentally obtained micrographs. Other common and effective approaches for microstructure modeling in- clude Markov random field (MRF) models whereby a Gibbs distribution is used to represent a microstructure and Gibbs sampling method is used for micrograph synthesis Geman and Geman (1984), Geman and Graffigne (1986). Parameters of the Markov random field model are typi- cally estimated from a small number of calibrating images and the micrographs are synthesized using Markov chain Monte Carlo (MCMC) methods. A large number of parameters is required to produce high resolution features with MRF models and estimating them with credible confidence from few calibrating images remains a challenge. Numerous variations on and enhancements to standard MRF have adapted it to specific situations which have included deterministic relaxation for image classification Berthod et al. (1996), Fan and Xia (2003), texture segmentation based on wavelets and hidden Markov tree models Choi and Baraniuk (2001). A Markov random field model for image segmentation that combines both color and texture features has also been pro- posed Kato and Pong (2000) as well as image imputing based on hierarchical Bayesian models to reconstruct image from partially observed multivariate data Dass and Nair (2003). Quantification of microstructural variance in the formalism of stochastic processes is developed in Niezgoda et al. (2011). Algorithms based on sequential testing and feature learning have also been de- veloped for face detection applications Fleuret and Geman (2001) while edge preserving image restoration with Huber-Markov random field models has been used to restore degraded images Pan and Reeves (2006), and procedures based on Gibbs reaction-diffusion have yielded a theory for learning generic prior models from a set of natural images Zhu and Mumford (1997). Maxi- mum a posteriori (MAP) estimation methods have also been used for image restoration Schultz and Stevenson (1994) sometimes in conjunction with adaptive parameter estimation Zhang et al. 83 (2005). An approach based on filters, random fields and maximum entropy (FRAME) Zhu et al. (1998) has also been very effective in capturing high resolution features in the image, but remains computationally expensive. Its equivalence with Julesz ensembles has also been investigated Wu et al. (2000). Finally, a parametric texture models based on an over complete complex wavelet transform was also used Portilla and Simoncelli (2000) to model micrographs based on the steer- able pyramid decomposition of texture (Simoncelli and Freeman (1995), Simoncelli et al. (2000)). In this paper, a similar texture model Portilla and Simoncelli (2000) is used with parameters ob- tained from a few high resolution calibrating images and low resolution image at the synthesis location. A particle filter is used to estimate the high resolution micrograph at a location given the low resolution micrograph at that location. In this work, a Bayesian framework is used, where the experimental low resolution micro- graph play the role of measurement used to update a prior model synthesized from a very small set of high resolution micrographs. A particle filter, namely a density-based Monte Carlo filter is used to implement the Bayesian formulation. In Tipireddy et al. (2009), Kalman filter is used for system identification based on multiple static and dynamic test data. Kalman filter methodology can not be used in our problem, because the state equation in the filtering problem is modeled as an implicit micrograph synthesis algorithm, where as the Kalman filter requires the state equation to be a linear equation. Monte Carlo based particle filters Doucet et al. (2000), Ghosh et al. (2008), Gordon et al. (1993), Manohar and Roy (2006), Tanizaki (1996) are used for system identification of dynamical system. In this work a density-based Monte Carlo filter proposed in Tanizaki (1996) is used. Synthesis of high resolution micrographs from low resolution images with parametric tex- ture model proposed in Portilla and Simoncelli (2000) along with the density-based Monte Carlo filter Tanizaki (1996) is demonstrated in our paper. The main contribution of the present work 84 is to compute the parameters of the texture model from high resolution calibrating images and low resolution images and combine particle filter with a texture synthesis model to characterize the high resolution micrograph. The paper is organized as follows. Motivation of the work is de- scribed with the problem statement in section “Problem statement”. In section “Parametric texture model”, the parametric model along with the synthesis algorithm is described. Section “Density- based Monte Carlo filter” illustrates the density-based Monte Carlo filter. Numerical experiments are provided in section “Numerical illustration”. 6.2 Problem statement Figure 6.1 shows a low resolution microscopic image of cast aluminum alloy W319. At locations A and B marked on the low resolution image, high resolution images are obtained from scanning electron microscopy. Figure 6.1: Low resolution micrograph Low resolution, low magnification digital image of cast aluminum alloy W319, for which high resolution images are available at points A and B 85 Figure 6.2 shows experimentally taken high resolution micrographs using scanning electron microscope at locations A and B. Figure 6.3 shows magnified low resolution micrographs avail- able at locations A, B and C. Our objective is to simulate high resolution micrographs at all the locations where the high resolution microscopic images are not available. The proposed high resolution micrograph synthesis methodology is demonstrated by generating a high resolution micrograph at location C. The procedure can be repeated to simulate the high resolution micro- graphs at other locations. (a) A (b) B Figure 6.2: High resolution digital image of cast aluminum alloy W319 at A and B Starting with an imagex 0 that represents the low-resolution micrograph at location C, our aim is to synthesize an imagex at the same location, representing the high-resolution micrograph, that is constrained, in some sense, by our knowledge of high resolution data at locations A and B and low resolution data at C. In order to complete the statement of the problem, we must specify the nature of the constraints (i.e. the specific functionals or features that describe our knowledge of the high resolution data). In what follows, a texture model will be used to provide a context in 86 (a) A (b) B (c) C Figure 6.3: Low resolution digital image of cast aluminum alloy W319 at A, B and C which features of a micrograph are quantified. This model will be constructed from available high resolution data and presumed to be valid everywhere else. Specifically, the texture model will be identified with the joint statistics of the subbands in a particular wavelet decomposition, namely, a steerable pyramid decomposition. Associated with this model is a mappingf that mapsx 0 into x such thatx has the joint statistics corresponding to the model. The texture model maps an initial micrograph x 0 into a high resolution micrograph x. Let x i 0 be a set of distinct initial micrographs, then the texture model maps each micrographx i 0 into a high resolution micrograph x i such that the parameters (joint statistics of the subbands) of each micrograph,x i are same as those of the texture model. However the visible microstructural features of each x i will be different from each other and also from those observed in the low resolution micrograph at that location. A Bayesian framework is employed using a particle filter, such that a high resolution micrographx is selected from all possible micrographsx i that is also consistent with the local microstructure. Based on maximum-entropy arguments, prior knowledge for the Bayesian framework takes the form of a Gaussian distribution with an interpolated low resolution micrograph as a mean and an assumed variance. The Likelihood is obtained from the 87 texture modelf, that maps low resolution imagex 0 to high resolution imagex. To complete the specification of the Bayesian problem, observations consisting of a low resolution micrograph available at location C are used to update the prior. A parametric texture model based on complex wavelet transform is used to represent the micrograph. High resolution micrograph is selected randomly from all possible micrographs satisfying a set of constraint functions, which will be discussed in section “Julesz conjecture and constraint functions”. The synthesis is performed by recursively projecting the initial guess of the micrograph onto the constraint surfaces using gradient projection methods Portilla and Simoncelli (2000). Density-based Monte Carlo filter Tanizaki (1996) along with the gradient projection algorithms are used to synthesize the high resolution micrograph at location C that is consistent with experimentally obtained low resolution micrograph at the same location. The problem is formulated as an estimation of a high resolution micrograph at location C using high resolution calibrating images available at locations A and B and low resolution image available at location C. The parametrization of the texture model is discussed in section “Parameter estimation sketch”. Texture parameters at location C are estimated from the parameters of images at locations A and B and those of coarse scale image at location C. The high resolution image at location C is simulated by recursively adjusting the pixel values of the iterates to satisfy a set of constraint functions and is updated using density-based Monte Carlo filter with coarse scale image at C as an observation. 88 6.3 Parametric texture model A texture is modeled as a real two dimensional random field,X(i;j) defined on a finite lattice, (i;j)2LZ 2 . HereX(i;j) takes gray level values in [0; 255]N. In this section we describe the steerable pyramid decomposition, which is a particular decomposition based on the complex wavelet transform of an image. A steerable pyramid acts as a filter and decomposes the image into multi-level, multi-oriented subbands. Joint statistics of the subbands, described in section “Pa- rameter estimation sketch” are used as the parameters of the texture model. Section “Micrograph synthesis procedure” describes the texture synthesis procedure. Julesz conjecture and constraint functions The Julesz conjecture (Julesz (1962), Zhu et al. (1998), Portilla and Simoncelli (2000)) states that there exists a set of functionsf k :R jLj 7!R;k = 1; ;N c g , called constraint functions Por- tilla and Simoncelli (2000) such that samples drawn from any two random fields are visually indistinguishable under some conditions, if expectations over this set of functions are equal. That is, E[ k (X)] =E[ k (Y )])XandY are perceptually indistinguishable; (6.1) where,f k ()g is the set of constraint functions acting on a random filed andE[] is mathematical expectation. To compute the expectations, practical ergodicity defined in Portilla and Simoncelli (2000) is assumed. Specifically, a homogeneous random field X has the property of practical ergodicity with respect to function : R jLj ! R, with tolerance ", and probability p, if and 89 only if the spatial average of over a realization drawn fromX is a good approximation to the expectation of with sufficiently high probability. More concretely, P x j(x(n;m))E[(X)]j<" p: (6.2) In the present case, where realizations refer to sample imagesx(i;j), the spatial average of the functional can be approximated as, (x(n;m)) = 1 L X (i;j)2L (x(bn +ic N ;bm +jc M )); (6.3) where,bc N means that the result is taken moduloN. With these definitions, our objective then is to draw samples from random fieldX satisfying statistical constraints of the form, E[ k (X)] =c k ;8k: (6.4) The valuesc k in these constraints are computed from calibrating imagex(i;j) and constitute the parameters of the texture model. The set of all samples satisfying equation 6.4 is referred to as the Julesz ensemble Zhu et al. (1998) given by T ;C =fx :E[ k (X)] =c k ;8kg: (6.5) 6.3.1 Steerable pyramid decomposition A parametric texture model based on over-complete complex wavelet transform Portilla and Si- moncelli (2000) is used in this work. The complex wavelet transform decomposes the calibrating 90 image into multi-scale, multi-oriented subbands. The joint statistics of the decomposed sub- bands are used as the parameters of the texture model. The decomposition is known as steerable pyramid when the basis functions are directional derivatives. The filters used in this decompo- sition are polar-separable in the Fourier domain. The traditional orthogonal separable wavelet decomposition suffers from aliasing whereas the steerable pyramid decomposition overcomes this limitation since the support of lowpass filter obeys the Nyquist sampling criterion Portilla and Simoncelli (2000). A disadvantage with the steerable pyramid decomposition is that the basis is over-complete. A system is said to be over-complete if the the number of basis functions used to represent a vector is larger than required. Translation and rotation invariance of the steerable pyramid is important to represent the oblique orientations properly. Filters used in the steerable pyramid belong to one of three categories, namely low-pass, high-pass or oriented band-pass filters which are characterized by their respective response functions, L(r) = 8 > > > > > > < > > > > > > : 2 cos 2 log 2 4r ; 4 <r< 2 2; r 4 0; r 2 ; (6.6) H(r) = 8 > > > > > > < > > > > > > : cos 2 log 2 2r ; 4 <r< 2 1; r 4 0; r 2 ; (6.7) and B k (r;) =H(r)G k (); k2 [0;K 1]; (6.8) 91 where the angular part ofB k is given by, G k () = 8 > > < > > : k cos k K K1 ; k K < 2 0; otherwise : (6.9) Here, r and are polar coordinates and k = 2 k1 (K1)! p K[2(K1)]! . L(r) and H(r) are low-pass and high-pass filters respectively andfB k (r;)g are oriented band-pass filters. The steerable pyramid decomposition is initiated by splitting an image into high-pass and low-pass subbands. The low-pass band is further decomposed into a set of orientation bands and the low-pass residual band. The low-pass residual band is sub-sampled by two and the decomposition in the next level is carried out by further decomposing the sub-sampled low-pass residual band into orientation bands and low-pass residual band. The procedure is continued till the image is decomposed into required number of pyramid levels. A system diagram for steerable pyramid is shown in fig 6.4. In our work four orientation bands and four pyramids with a total of 18 subbands are used (16 orientation bands, one high-pass and one low pass residual band). Figure 6.5 shows the real parts of four level and four orientation steerable pyramid subbands and low-pass residual band. High-pass subband is not shown in figure 6.5. 6.3.2 Parameter estimation sketch Outline for estimating parameters of the texture model from a sample image is described in the following steps: (i). Decompose the high resolution sample image x into high-pass subband, h 0 , orientation bands,fb s;k ;s = 1; ; 4;k = 1; ; 4g and low-pass residual bands,fl 0 ;l 1 ;l 2 ;l 3 ;l 4 g. 92 Figure 6.4: Steerable pyramid diagram Block diagram for the steerable pyramid Portilla and Simoncelli (2000), Simoncelli and Freeman (1995). Input image is initially split into high-pass and low-pass bands. Low-pass band is further decomposed into a set of oriented subbands and low-pass residual band. The recursive decomposition is achieved by inserting the diagram in dotted lines into the solid circle. Figure 6.5: Real parts of subbands Real parts of steerable pyramid orientation bands(4 orientations and 4 scales) and low-pass subband of image at location A. High-Pass subband is not shown here. 93 (ii). Compute a set of joint statistics on the decomposed subbands. Following are the statis- tics Portilla and Simoncelli (2000) that are considered as parameters of the texture model: Marginal statistics are the statistics of gray-level values of the pixels in the image. – Mean, 1 (x) =Efx(i;j)g. – Variance, 2 (x) =Ef(x(i;j) 1 (x)) 2 g – Skewness,(x) = 3 (x) ( 2 (x)) 1:5 – Kurtosis,(x) = 4 (x) ( 2 (x)) 2 – Range = [min(x)); max(x)] where, n (x) = Ef(x(i;j) 1 (x)) n g;8n > 1. In the parametric texture model, mean 1 (x), variance 2 (x), skewness(x), kurtosis(x), minimum and maximum of the sample image, variance of the high-pass band 2 (h 0 ), skewnessf(l 0 );(l 1 );(l 2 );(l 3 );(l 4 )g and kurtosis,f(l 0 );(l 1 );(l 2 );(l 3 );(l 4 )g of partially reconstructed low-pass im- ages at each scale,fl 0 ;l 1 ;l 2 ;l 3 ;l 4 g are considered as a set of parameters from marginal statistics. Raw coefficient correlations, are defined as the central samples of auto-correlation of the partially reconstructed low-pass images,A k (n;m) =Efl k (i;j)l k (bi+nc N ;bj + mc M )g. In the texture model, auto correlations of the partially reconstructed low-pass images,fl 0 ;l 1 ;l 2 ;l 3 ;l 4 g at each scale, are used as a set of parameters. These statistics characterize the salient spatial frequencies and the regularity of the texture Portilla and Simoncelli (2000). Coefficient magnitude statistics are defined as the central samples of auto-correlation of magnitude of each subband, cross correlation of magnitude of each subband with 94 those of the other orientations at the same scale, cross correlation of subband mag- nitudes with those of all orientations at a coarser scale. Cross-correlation of any two subbandsb s (i;j) andb t (i;j) is given asC s;t (n;m) = Efb s (i;j)b t (bi +nc N ;bj + mc M )g. These represent the structures such as edges, bars and corners in the image. Cross-scale phase statistics are cross-correlation of the real part of the coefficients with both the real and imaginary part of the phase doubled coefficients at all orien- tations at the next coarser scale. The rate at which the phase changes for fine scale coefficients is twice the rate change for coarse scale coefficients Portilla and Simon- celli (2000). To compensate for the difference in rate of phase change, the cross correlations are computed between the phase doubled coarse scale coefficients and fine scale coefficients. These statistics distinguish edges from lines. All the joint statistics described in the above list are used as the parameters of the texture model. Practical ergodicity described in section “Julesz conjecture and constraint func- tions” is used to compute expectationsE[] in the above expressions . (iii). Once a set of parameters is computed from high resolution micrograph, a micrograph with that set of parameters can be constructed using the synthesis procedure described in the next section. 6.3.3 Micrograph synthesis procedure The synthesis procedure proposed in Portilla and Simoncelli (2000) to generate an image with desired parameters is described in the following steps. The parameters of the texture model are computed from calibrating images as joint statistics of the decomposed subbands. 95 (i). The synthesis procedure is initiated by choosing as initial guess an image with mean and variance equal to those of the target image. (ii). The initial image is decomposed into multi-scale, multi-oriented subbands and joint statis- tics of the desired image are imposed recursively using gradient projection method Portilla and Simoncelli (2000). (iii). Starting from coarse scale to fine scale, statistical constrains are imposed on low-pass and orientation bands while simultaneously reconstructing the low-pass image. (iv). Auto-correlation, skewness and kurtosis of the reconstructed low-pass image at the fine scale are adjusted and resulting low-pass image is added to variance adjusted high-pass image to obtain the fine scale image. (v). The marginal statistics are adjusted on the pixels of the reconstructed image in step-4. (vi). The procedure is repeated from step-2 to step-5 until the image with desired statistics is generated. (vii). This algorithm synthesizes an image that is statistically equivalent to the sample images from which the parameters are estimated. Fig 6.6 shows one iteration step of synthesis algorithm. This algorithm can be used to simulate many statistically equivalent micrographs for a given set of parameters. But the simulated micro- graphs will not be visually similar to the calibrating micrograph from which the parameters are computed. When a low resolution microscopic image from an experiment is available, data as- similation methods can be used to synthesize a high resolution micrograph that is consistent with the low resolution image. We use the Matlab tools provided at Simoncelli (2009) for steerable 96 pyramid decomposition and micrograph synthesis. In the next section a particle filter known as the density-based Monte Carlo filter Tanizaki (1996) along with micrograph synthesis procedure is used to simulate high resolution micrographs that are consistent with the corresponding low resolution one. Figure 6.6: Block diagram of the recursive micrograph synthesis procedure 6.4 Density-based Monte Carlo filter The density-based Monte Carlo filter is a nonlinear, non-Gaussian filter proposed in Tanizaki (1996) for state estimation of a dynamical system. In the current problem, an image is synthesized by an extension of the procedure described in section “Micrograph synthesis procedure”. In the synthesis procedure, an initial guess of the image is assumed and the desired joint statistics are imposed on the initial guess to obtain a new image. Letf be an implicit algorithm that implements steps 2-5 in the synthesis procedure shown in section “Micrograph synthesis procedure”. The parameters of the texture model are repeatedly imposed on the initial image x 0 such that the final image has the desired parameters. That is the implicit function f acts on x 0 to generate x 1 , on x 1 to generate x 2 and in general acts on x k1 to generate x k . Action of f on x k1 to generate x k , is shown in figure 6.6. We treat this as one function evaluation in the synthesis algorithm. Let x 0 ;x 1 ; ;x k denote the images, modeled as random fields, generated at each function evaluation of the synthesis algorithm. These are construed as successive states of a 97 dynamical system where an implicit algorithm takes a dynamical state (in the present case, the state is a random field) x k1 to x k where x k depends only on x k1 and the parameters of the texture model. Here, function evaluation steps, k = 0; ;K are treated analogous to time steps in the dynamical system. Coarse scale images available at synthesis locations are treated as measurements and are used in the measurement equation of the particle filter algorithm to update the estimated image at each synthesis step. The state and measurement equations of the dynamical system can thus be written as x k = f(x k1 ) + k ; state equation (6.10) y k = g(x k ) +" k ; measurement equation (6.11) where,x k is the synthesized image atk th step in the recursive synthesis process,f is an implicit algorithm relatingx k tox k1 as shown in figure 6.6 andg is a coarsening operator relating state x k and observationy k . Here, coarsening operator is an averaging operator that averages 8 8 pixels in fine scale image to produce one pixel in the coarse scale image. Errors in the synthesis process and experimental measurements are modeled as additive white Gaussian noise processes k and " k respectively. The covariance of the noises are assumed to be known and are indeed parameters of the filter. Let X k =fx 0 ; ;x k g be a vector of states and Y k =fy 0 ; ;y k g be a vector of measurements. In our work, measurements are assumed to consist of a coarse scale image available at location C and does not change with each synthesis step. Hence, we assume,Y k =fy 0 ;y 1 =y 0 ; ;y k =y 0 g. The joint probability density of all the states and the measurements is given by P (X k ;Y k ) =P (Y k jX k )P (X k ); (6.12) 98 where,P (Y k jX k ) andP (X k ) can be written as, P (Y k jX k ) = k Y s=1 P (y s jx s ) (6.13) and P (X k ) =P (x 0 ) k Y s=1 P (x s jx s1 ): (6.14) The conditional densities, P (y s jx s ) and P (x s jx s1 ), are obtained from the measurement and state equations respectively. From equations 6.13 and 6.14, the filter density P (x k jY k ) can be written as, P (x k jY k ) = R P (Y k jX k )P (X k )dX k1 R P (Y k jX k )P (X k )dX k : (6.15) The estimated micrograph, ^ x k at step k, can be obtained by computing expectation of x k with respect to the filtering density,P (x k jY k ), i.e., ^ x k = R x k P (Y k jX k )P (X k )dX k R P (Y k jX k )P (X k )dX k : (6.16) Evaluating the integrals in equation 6.16 using Monte Carlo sampling results in the filter estimate ^ x k at stepk. Generatingn random samples ofx 0 and s fors =f1; ;kg, a set of random drawsX i;s can be obtained from the state equation as x i;s =f(x i;s1 ; i;s ); i = 1; ;n &s = 1; ;k: (6.17) 99 The filtering estimate ^ x k in equation (6.16) is approximated by drawing samples of x k using Monte Carlo sampling according to, ^ x k = P n i=1 x i;k P (Y k jX i;k ) P n i=1 P (Y k jX i;k ) : (6.18) Using equation 6.12, ^ x k = P n i=1 x i;k Q k s=1 P (y s jx i;s ) P n i=1 Q k s=1 P (y s jx i;s ) ; (6.19) where, (y s jx i;s ) N ((g s (x i;s )); cov(" s )) denotes the likelihood and ^ x k is the filter estimate. Here,X i;k is a collection of random draws, defined asX i;k =fx i;0 ; ;x i;k g. Equation (6.18) can be written as ^ x k = n X i=1 x i;k w i;k ; (6.20) where the weight functionw i;k is defined as w i;k = P (y k j(x i;k )w i;k1 P n j=1 P (y k j(x j;k )w j;k ; (6.21) such that P n i=1 w i;k = 1. For the initial step, whenk = 0, equal weights are assigned for all the samplesfx i ; 0g n i=1 , i.e.,w i;0 = 1 n . Writinge k = (x k ^ x k ), the variance ofe k can be written as k = R e k e T k (P (Y k jX k )) 2 P (X k )dX k ( R (P (Y k jX k ))P (X k )dX k ) 2 (6.22) and the sample variance can be written as k = P n i=1 e k e T k (P (Y k jX k )) 2 P (X k ) ( P n i=1 (P (Y k jX k ))P (X k )) 2 (6.23) 100 k = n X i=1 (x i;k ^ x k )(x i;k ^ x k ) T w 2 i;k : (6.24) The following limiting conditions asn goes to infinity can be verified Tanizaki (1996) ^ x k !E(x k jY k ) almost surely, p n(x k E(x k jY k ))!N(0; k ) in distribution , n k ! k almost surely. Following the above constructions, the computational procedure of the density-based Monte Carlo filter is as follows: (i). The random draws of the initial statex 0 , represented asx i;0 are taken from the initial density P (x 0 ), (ii). Givenx i;k1 and i;k1 , estimates ofx i;k are obtained from the state equation 6.10, (iii). Given the initial weightw i;0 , weight functionsw i;k fori = 1; ;n andk = 1; ;K are obtained from equation 6.21, (iv). Finally, the filtered state ^ x k is evaluated from equation 6.20 while the variance of the error in the estimate is computed from equation 6.24. Fig 6.7 shows the block diagram of the micrograph synthesis using parametric texture model and density-based Monte Carlo filter. 6.5 Numerical illustration The high resolution micrograph synthesis procedure is demonstrated by synthesizing a high res- olution micrograph at location C starting with low resolution micrograph at that location and 101 Figure 6.7: Block diagram of synthesis procedure with particle filter. Block diagram of the recur- sive micrograph synthesis procedure along with the density-based Monte Carlo filter. high resolution calibrating micrographs available at locations A and B as shown in figures 6.1 and 6.2. The procedure can be repeated to construct high resolution micrographs at any other location. Two cases are considered to construct high resolution micrographs. In the first case (figures 6.8, 6.10 and 6.12), a parametric texture model with the synthesis procedure described in section “Micrograph synthesis procedure” is used to construct a high resolution micrograph, whereas in the second case (figures 6.9, 6.11 and 6.13), the parametric texture model along with the density-based Monte Carlo filter described in section “Density-based Monte Carlo filter” is used to estimate a high resolution micrograph. (a) High resolution A (b) micrograph-1 (c) micrograph-2 (d) micrograph-3 Figure 6.8: Micrograph at A synthesized without particle filter. Micrographs from left to right show respectively, the experimental high resolution micrograph, synthesized micrographs without particle filter which are statistically equivalent with high resolution micrograph at location A. 102 (a) High resolution A (b) Low resolution A (c) Estimated A Figure 6.9: Micrograph at A synthesized with particle filter. From 6.9(a) to 6.9(c), experimental high resolution micrograph, experimental low resolution micrograph and high resolution micro- graph synthesized with texture model and particle filter at location A. (a) High resolution B (b) micrograph-1 (c) micrograph-2 (d) micrograph-3 Figure 6.10: Micrograph at B synthesized without particle filter. Micrographs from left to right show respectively, the experimental high resolution micrograph, synthesized micrographs without particle filter which are statistically equivalent with high resolution micrograph at location B. When the particle filter is not used, the synthesis procedure does not make use of the low res- olution micrographs available at a given location. It uses only the parameters of the texture model obtained from the calibrating micrographs. The synthesis process with distinct initial guesses re- sults in statistically equivalent high resolution micrographs, that have visually distinct microstruc- tural features. To validate the synthesis procedure, high resolution micrographs are synthesized at locations A and B where experimental high resolution micrographs are available. Figure 6.8(a) 103 (a) High resolution B (b) Low resolution B (c) Estimated B Figure 6.11: Micrograph at B synthesized with particle filter. From6.11(a) to 6.11(c), experi- mental high resolution micrograph, experimental low resolution micrograph and high resolution micrograph synthesized with texture model and particle filter at location B. (a) High resolution C (b) micrograph-1 (c) micrograph-2 (d) micrograph-3 Figure 6.12: Micrograph at C synthesized without particle filter. Micrographs from left to right show respectively, the experimental high resolution micrograph, synthesized micrographs without particle filter which are statistically equivalent with high resolution micrograph at location C. High resolution micrograph at location C is used only for comparison and is not used in the synthesis process. shows a high resolution micrograph obtained experimentally at location A and figures 6.8(b) – 6.8(d) show synthesized micrographs whose joint statistics are the same as those of 6.8(a) but with visibly distinct microstructural features. Similarly figure 6.10(a) shows an experimental high resolution micrograph at location B and figures 6.10(b) – 6.10(d) show corresponding syn- thesized micrographs. Parameters of the texture model for the synthesis of these micrographs at locations A and B are obtained using high resolution micrographs available at these locations. 104 (a) High resolution C (b) Low resolution C (c) Estimated C Figure 6.13: Micrograph at C synthesized with particle filter. From 6.13(a) to 6.13(c), experi- mental high resolution micrograph, experimental low resolution micrograph and high resolution micrograph synthesized with texture model and particle filter at location C. High resolution mi- crograph at location C is used only for comparison and is not used in the synthesis process. In case of location C, on the other hand, we assume that an experimental high resolution micro- graph is not available and we aim to synthesize one. The parameters of the texture model for C are obtained from a low resolution micrograph at C and high resolution micrographs available at locations A and B. Pixel statistics and low-pass band statistics which account for total gray scale levels are computed from low resolution micrograph at C and the rest of the parameters which determine high resolution feature are computed as averages of the parameters from A and B. Figure 6.12(a) shows experimental high resolution micrograph and figures 6.12(b) – 6.12(d) show corresponding synthesized micrographs at location C. Here high resolution micrograph at location C (figure 6.12(a)) is used only for comparison and is not used in the synthesis process. In the second case, high resolution micrographs of the microstructure are estimated in the Bayesian context using the density-based Monte Carlo filter (section “Density-based Monte Carlo filter”) where the texture synthesis process (section “Micrograph synthesis procedure”) is treated as a dynamical system and the experimental low resolution micrograph at a given location is 105 treated as measurement used to update the dynamical state. It is reiterated that the high resolu- tion micrograph at C is synthesized using the low resolution micrograph at C and high resolution micrographs at A and B together with a texture model and a particle filter. The parameters of the texture model for A, B and C are obtained similarly to before. For estimating high resolution micrographs at locations A, B and C, experimental low resolution micrographs at corresponding locations are used as measurements in the filtering algorithm. Low resolution micrographs are also used as the statistical average of the initial probability density function in the algorithm. The algorithm convergence is slowed for other choices of initial average. Measurement and process noises are assumed to be zero mean Gaussian random fields with 0:01 coefficient of variation. In the particle filter algorithm, 25 Monte Carlo samples are used to estimate the high resolu- tion micrographs. Figure 6.9(a) and figure 6.9(b) show high and low resolution micrographs obtained experimentally at location A, figure 6.9(c) is the estimated high resolution micrograph at location A using density-based Monte Carlo filter along with the texture synthesis procedure. Figures 6.11(a)– 6.11(c) show the experimental high resolution, experimental low resolution, es- timated high resolution micrographs at location B. At location C, figure 6.13(a) and figure 6.13(b) are experimental high and low resolution micrographs, figure 6.13(c) is the estimated high reso- lution micrograph. It should be noted that the high resolution features of micrograph at location C (figure 6.13(c)), are recovered using information from experimental low resolution micrograph at C and experimental high resolution micrographs at A and B. It is again worth comparing figure 6.12 showing the synthesis of high resolution micrograph using only the parametric texture model with figure 6.13 showing high resolution micrograph synthesized using texture model along with the particle filter. It can be clearly observed that the estimated high resolution micrograph with particle filter (figure 6.13(c)) provides information 106 about the local microstructure that cannot be gleaned from previous approaches. Resolving this level of detail is paramount in physical problems where bifurcation is initiated by spatially local events. 6.6 Appendix Here, we provide the pseudo-code for high resolution micrograph synthesis using the parametric texture model and density based Monte-Carlo filter. A few high resolution micrographs are de- composed using steerable pyramid method and the parameters of the texture model are obtained as joint statistics of the decomposed subbands described in “Parameter estimation sketch”. Algorithm 3 Recursive micrograph syntheis procedure 1. Compuete parametersp from sample micrographs 2. Assume initial guessx 0 3. Compute parameters ~ p fromx 0 4. while ~ p6=p do 5. Decomposex k into subbands,h 0 ,fb s;k ;s = 1; ; 4;k = 1; ; 4g andfl 0 ;l 1 ;l 2 ;l 3 ;l 4 g 6. Compute the parameters ~ p from above subbands 7. From coarse to fine scale, impose desired statistics and reconstructl 0 8. Impose auto-correlation onl 0 9. Impose skewness and kurtosis onl 0 10. Impose variance onh 0 11. Construct finescale micrographx k+1 =l 0 +h 0 12. Impose pixel level statistics onx k+1 13. end while. 107 Algorithm 4 Density based Monte-Carlo filter 1. Draw initial samplesx i;0 from densityP (x 0 ) 2. fork = 1 toK do 3. Computex i;k =f(x i;k1 + i;k1 ), wheref is micrograph synthesis algorinthm 3 4. Compute wightsw i;k from equation 6.21 5. end for 6. Estimate ^ x from equation 6.20 and variance from equation 6.24. 108 Bibliography Adams, B. M., Dalbey, K. R., Eldred, M. S., Gay, D. M., Swiler, L. P., Bohnhoff, W. J., Eddy, J. P., Haskell, K., and Hough, P. D. 2010. DAKOTA, A Multilevel Parallel Object-Oriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis. Sandia National Laboratories, tech. rep.sand2010-2183 edition. Babuˇ ska, I., Tempone, R., and Zouraris, G. E. 2004. Galerkin finite element approximations of stochastic elliptic partial differential equations. SIAM J. Numer. Math., 42:800–825. Berthod, M., Kato, Z., Yu, S., and Zerubia, J. 1996. Bayesian image classification using markov random fields. Image and Vision Computing, 14:285–295. Betchen, L., Straatman, A. G., and Thompson, B. E. 2006. A nonequilibrium finite-volume model for conjugate fluid/porous/solid domains. Numerical Heat Transfer, Part A, 49:543–565. Brooks, A. N. and Hughes, T. J. 1982. Streamline upwind/petrov-galerkin formulations for con- vection dominated flows with particular emphasis on the incompressible navier-stokes equations. Computer Methods in Applied Mechanics and Engineering, 32:199–259. C. Gear, J. H., Kevrekidis, P., Kevrekidis, I., Runborg, O., and Theodoropoulos, C. 2003. Equation-free, coarse-grained multiscale computation: Enabling mocroscopic simulators to per- form system-level analysis. Commun. Math. Sci., 1:715–762. Chandesris, M., Serre, G., and Sagaut, P. 2006. A macroscopic turbulence model for flow in porous media suited for channel, pipe and rod bundle flows. International Journal of Heat and Mass Transfer, 49:27392750. Choi, H. and Baraniuk, R. G. 2001. Multiscale image segmentation using wavelet-domain hid- den markov models. IEEE Transactions on Image Processing, 10:1309–1321. Dass, S. C. and Nair, V . N. 2003. Edge detection, spatial smoothing, and image reconstruction with partially observed multivariate data. Journal of the American Statistical Association, 98:77– 89. Doostan, A., Ghanem, R., and Red-Horse, J. 2007. Stochastic model reduction for chaos repre- sentations. Computer Methods in Applied Mechanics and Engineering, 196:3951–3966. Doostan, A. and Iaccarino, G. 2009. A least-squares approximation of partial differential equa- tions with high-dimensional random inputs. Journal of Computational Physics, 228:4332–4345. 109 Doostan, A. and Owhadi, H. 2011. A non-adapted sparse approximation of pdes with stochastic inputs. Journal of Computational Physics, 230:3015–3034. Doucet, A., Godsill, S., and Andrieu, C. 2000. On sequential monte carlo sampling methods for bayesian filtering. Statistics and computing, 10:197–208. Elman, H. C., Miller, C. W., Phipps, E. T., and Tuminaro, R. S. 2011. Assessment of collocation and galerkin approaches to linear diffusion equations with random data. International Journal for Uncertainty Quantifications, 1:19–34. Fan, G. and Xia, X. G. 2003. Wavelet-based texture analysis and synthesis using hidden markov models. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 50:106–120. Fleuret, F. and Geman, D. 2001. Coarse-to-fine face detection. International Journal of Com- puter Vision, 41:85–107. Geman, S. and Geman, D. 1984. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Pat. Anal. Mach. Intell, 6:721–741. Geman, S. and Graffigne, C. 1986. Markov random field image models and their applications to computer vision. In Proceedings of International Congress of mathematicians, pages 1496– 1517, Berkeley, California, USA. Ghanem, R. 1999. Ingredients for a general purpose stochastic finite elements implementation. Computer Methods in Applied Mechanics and Engineering, 168:19–34. Ghanem, R. and Spanos, P. 1991. Stochastic Finite Elements: A Spectral Approach. Springer- Verlag. Ghanem, R. G. and Doostan, A. 2006. On the construction and analysis of stochastic models: characterization and propagation of the errors associated with limited data. Journal of Compu- tational Physics, 213:63–81. Ghosh, S., Valiveti, D., Harris, S. J., and Boileau, J. 2006. A domain partitioning based pre- processor for multi-scale modelling of cast aluminium alloys. Modelling and simulation in materials science and engineering, 14:1363–1396. Ghosh, S. J., Manohar, C. S., and Roy, D. 2008. A sequential importance sampling filter with a new proposal distribution for state and parameter estimation of nonlinear dynamical systems. Proc. R. Soc. A, 464:25–47. Gordon, N. J., Salmond, D. J., and Smith, A. F. M. 1993. Novel approach to nonlinear/non- gaussian bayesian state estimation. IEEE Trans. Radar Signal Process., 140:107–113. Groeber, M., Ghosh, S., Uchic, M. D., and Dimiduk, D. M. 2008a. A framework for automated analysis and simulation of 3d polycrystalline microstructures.: Part 1: Statistical characteriza- tion. Acta Materialia, 56:12571273. 110 Groeber, M., Ghosh, S., Uchic, M. D., and Dimiduk, D. M. 2008b. A framework for auto- mated analysis and simulation of 3d polycrystalline microstructures. part 2: Synthetic structure generation. Acta Materialia, 56:12741287. Holden, H., Oksendal, B., Uboe, J., and Zhang, T. 1996. Stochastic Partial Differential Equa- tions: a modeling, white noise functional approach. Birkhauser Boston Inc., Cambridge, MA, USA. Janson, S. 1999. Gaussian Hilbert Spaces. Cambridge University Press. Julesz, B. 1962. Visual pattern discrimination. IRE Trans. Info Theory, IT-8:84–92. Kato, Z. and Pong, T. C. 2000. A markov random field image segmentation model for color textured images. Image and Vision Computing, 24:1103–1114. Lehmann, T. M., G¨ onner, C., and K.Spitzer 1999. Survey: Interpolation methods in medical image processing. IEEE Transactions on Medical Imaging, 18:1049–1075. Manohar, C. S. and Roy, D. 2006. Monte carlo filters for identification of nonlinear structural dynamical systems. Sadhana, 31(4):399–427. Nakayama, A., Kuwahara, F., Umemoto, T., and Hayashi, T. 2002. Heat and fluid flow within an anisotropic porous medium. Transactions of ASME, 124:746–753. Niezgoda, S. R., Yabansu, Y . C., and Kalidindi, S. R. 2011. Understanding and visualizing mi- crostructure and microstructure variance as a stochastic process. Acta Materialia, 59:63876400. Nobile, F., Tempone, R., and Webster, C. G. 2008. A sparse grid stochastic collocation method for partial differential equations with random input data. SIAM Journal of Numerical Analysis, 46:2309–2345. Pan, R. and Reeves, S. J. 2006. Efficient huber-markov edge-preserving image restoration. IEEE Transactions on Image Processing, 15:3728–3735. Portilla, J. and Simoncelli, E. 2000. A parametric texture model based on joint statistics of complex wavelet coefficients. International Journal of Computer Vision, 40(1):49–70. Powell, C. and Elman, H. 2009. Block-diagonal preconditioning for spectral stochastic finite element systems. IMA Journal of Numerical Analysis, 29:350–375. Rollett, A. D., Lee, S. B., Campman, R., and Rohrer, G. 2007. Three-dimensional characteriza- tion of microstructure by electron back-scatter diffraction. Annual Review of Materials Research, 37:627–658. Rosseel, E. and Vandewalle, S. 2010. Iterative solvers for the stochastic finite element method. SIAM J. Sci. Comput, 32:372–397. Saad, Y . 1996. Iterative Methods for Sparse Linear Systems. SIAM. Schmidt, R. C., Belcourt, K., Clarno, K. T., Hooper, R., Humphries, L. L., Lorber, A. L., Pryor, R. J., and Spotz, W. F. 2010. Foundational Development of an Advanced Nuclear Reactor Integrated Safety Code. Sandia National Labs, sand2010-0878 edition. 111 Schultz, R. and Stevenson, R. L. 1994. A bayesian approach to image expansion for improved definition. IEEE Transactions on Image Processing, 3:233–242. Schwab, C. and Gittelson, C. 2011. Sparse tensor discretizations of high-dimensional parametric and stochastic pdes. Acta Numerica, 20:291–467. Shadid, J., Salinger, A., Pawlowski, R., Lin, P., Hennigan, G., Tuminaro, R., and Lehoucq, R. 2006. Large-scale stabilized fe computational analysis of nonlinear steady-state trans- port/reaction systems. Comput. Methods Appl. Mech. Engrg., 195:1846–1871. Simoncelli, E. P. 2009. matlabpyrtools code. Simoncelli, E. P. and Freeman, W. 1995. The steerable pyramid: A flexible architecture for multi-scale derivative computation. In In Second Int’l Conf. on Image Proc., volume III, pages 444–447, Washington, DC,. IEEE Sig. Proc. Society. Simoncelli, E. P., Freeman, W., Adelson, E., and Heeger, D. 2000. Shiftable multi-scale trans- forms. IEEE Trans Information Theory, 38(2):587–607. Singh, H., Gokhale, A. M., Mao, Y ., and Spowart, J. E. 2006. Computer simulations of realistic microstructures of discontinuously reinforced aluminum alloy (dra) composites. Acta Materi- alia, 54:21312143. Smolyak, S. 1963. Quadrature and interpolation formulas for tensor products of certain classes of functions. Dokl. Akad. Nauk SSSR, 148:1042–1043. Tanizaki, H. 1996. Nonlinear filters: estimation and applications. Springer Verlag, Berlin, second edtion edition. Tewari, A., Gokhale, A. M., Spowart, J. E., and Miracle, D. B. 2004. Quantitative charac- terization of spatial clustering in three-dimensional microstructures using two-point correlation functions. Acta Materialia, 52:307319. Tezduyar, T. E. 1992. Stabilizedfinite element formulations for incompressibleflow computation. Advances in Applied Mechanics, 28:1–44. Tipireddy, R., Nasrellah, H. A., and Manohar, C. S. 2009. A kalman filter based strategy for lin- ear structural system identification based on multiple static and dynamic test data. Probabilistic Engineering Mechanics, 24:60–74. Ullmann, E. 2010. A kronecker product preconditioner for stochastic galerkin finite element discretizations. SIAM J. Sci. Comput., 32(2):923–946. Valiveti, D. M. and Ghosh, S. 2007. Morphology based domain partitioning of multi-phase materials: A preprocessor for multi-scale modelling. Int. J. Numer. Meth. Engng, 69:1717–1754. Weiner, N. 1963. The homogeneous chaos. American Journal of Mathematics, 60:897–936. Wu, Y . N., Zhu, S. C., and Liu, X. 2000. Equivalence of julesz ensembles and FRAME models. International Journal of Computer Vision, 38:247–265. 112 Xiu, D. and Karniadakis, G. 2002. The wiener-askey polynomial chaos for stochastic differential equations. SIAM J. Sci. Comput., 24:619–644. Zhang, X., Lam, K. M., and Shen, L. 2005. Image magnification based on adaptive mrf model parameter estimation. In Proceedings of 2005 International Symposium on Intelligent Signal Proceeding and Communication Systems, pages 653–656, Hong Kong. Zhu, S., Wu, Y . N., and Mumford, D. 1998. Filters, random fields and maximum entropy (FRAME):towards the unified theory for texture modeling. International Journal of Computer Vision, 27(2):107–126. Zhu, S. C. and Mumford, D. 1997. Prior learning and gibbs reaction-diffusion. IEEE Transac- tions on Pattern Analysis and Machine Intelligence, 19:1236–1250. 113

Abstract (if available)

Abstract This dissertation focuses on facilitating the analysis of probabilistic models for physical systems. To that end, novel contributions are made to various aspects of the problem, namely, 1) development of efficient algorithms for solving stochastic system of equations governing the physical problem, 2) stochastic basis adaptation methods to compute the solution in reduced stochastic dimensional space and 3) stochastic upscaling methods to find coarse-scale models from the fine-scale stochastic solutions. In particular, algorithms are developed for stochastic systems that are governed by partial differential equations (PDEs) with random coefficients. Polynomial chaos-based stochastic Galerkin and stochastic collocation methods are employed for solving these equations. Solvers and preconditioners based on Gauss-Seidel and Jacobi algorithms are explored for solving system of linear equations arising from stochastic Galerkin discretization of PDEs with random input data. Gauss-Seidel and Jacobi algorithms are formulated such that the existing software is leveraged in the computational effort. These algorithms are also used to develop preconditioners to Krylov iterative methods. These solvers and preconditioners are tested by solving a steady state diffusion equation and a steady state advection-diffusion equation. Upon discretization, the former PDE results in a symmetric positive definite matrix on left-hand-side, whereas the latter results in a non-symmetric positive definite matrix. The stochastic systems face significant computational challenge due the curse of dimensionality as the solution often lives in very high dimensional space. This challenge is addressed in the present work by recognizing the low dimensional structure of many quantities of interest (QoI) even in problems that have been embedded, via parameterization, in very high-dimensional settings. A new method for the characterization of subspaces associated with low-dimensional QoI is presented here. The probability density function of these QoI is found to be concentrated around one-dimensional subspaces for which projection operators are developed. This approach builds on the properties of Gaussian Hilbert spaces and associated tensor product spaces. ❧ For many physical problems, the solution lives in multiple scales, and it is important to capture the physics at all scales. To address this issue, a stochastic upscaling methodology is developed in which the above developed algorithms and basis adaptation methods are used. In particular upscaling methodology is demonstrated by developing a coarse scale stochastic porous medium model that replaces a fine-scale which consists of flow past fixed solid inclusions. The inclusions have stochastic spatially varying thermal conductivities and generate heat that is transported by the fluid. The permeability and conductivity of the effective porous medium are constructed as statistically dependent stochastic processes that are both explicitly dependent on the fine scale random conductivity. ❧ Another contribution of this thesis is development of a probabilistic framework for synthesizing high resolution micrographs from low resolution ones using a parametric texture model and a particle filter. Information contained in high resolution micrographs is relevant to the accurate prediction of microstructural behavior and the nucleation of instabilities. As these micrographs may be tedious and uneconomical to obtain over an extended spatial domain, A statistical approach is proposed for interpolating fine details over a whole computational domain starting with a low resolution prior and high resolution micrographs available only at a few spatial locations. As a first step, a small set of high resolution micrographs are decomposed into a set of multi-scale and multi-orientation subbands using a complex wavelet transform. Parameters of a texture model are computed as the joint statistics of the decomposed subbands. The synthesis algorithm then generates random micrographs satisfying the parameters of the texture model by recursively updating the gray level values of the pixels in the input micrograph. A density-based Monte Carlo filter is used at each step of the recursion to update the generated micrograph, using a low resolution micrograph at that location as a measurement. The process is continued until the synthesized micrograph has the same statistics as those from the high resolution micrographs. The proposed method combines a texture synthesis procedure with a particle filter and produces good quality high resolution micrographs.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

PDF

Comprehensive uncertainty quantification in composites manufacturing processes

PDF

Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

PDF

Model, identification & analysis of complex stochastic systems: applications in stochastic partial differential equations and multiscale mechanics

PDF

Stochastic and multiscale models for urban and natural ecology

PDF

Risk assessment, intrinsic interpolation and computationally efficient models for systems under uncertainty

PDF

Efficient stochastic simulations of hydrogeological systems: from model complexity to data assimilation

PDF

A polynomial chaos formalism for uncertainty budget assessment

PDF

Inverse modeling and uncertainty quantification of nonlinear flow in porous media models

PDF

Hybrid physics-based and data-driven computational approaches for multi-scale hydrologic systems under heterogeneity

PDF

Machine-learning approaches for modeling of complex materials and media

PDF

Modeling of multiscale continuum-atomistic systems using homogenization theory

PDF

Probabilistic data-driven predictive models for energy applications

PDF

Second order in time stochastic evolution equations and Wiener chaos approach

PDF

Stochastic models: simulation and heavy traffic analysis

PDF

Stochastic oilfield optimization under uncertainty in future development plans

PDF

Computational validation of stochastic programming models and applications

PDF

Computational stochastic programming with stochastic decomposition

PDF

Algorithms and landscape analysis for generative and adversarial learning

PDF

Studies into computational intelligence approaches for the identification of complex nonlinear systems

Asset Metadata

Creator Tipireddy, Ramakrishna (author)

Core Title Algorithms for stochastic Galerkin projections: solvers, basis adaptation and multiscale modeling and reduction

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering

Publication Date 08/02/2013

Defense Date 06/20/2013

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag curse of dimensionality,micrograph synthesis,model reduction,OAI-PMH Harvest,polynomial chaos,stochastic analysis,stochastic upscaling,uncertainty quantifi cation

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ghanem, Roger (committee chair), Lucas, Robert (committee member), Masri, Sami (committee member), Nakano, Aiichiro (committee member)

Creator Email rktreddy@gmail.com,tipiredd@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-311834

Unique identifier UC11294437

Identifier etd-TipireddyR-1935.pdf (filename),usctheses-c3-311834 (legacy record id)

Legacy Identifier etd-TipireddyR-1935.pdf

Dmrecord 311834

Document Type Dissertation

Format application/pdf (imt)

Rights Tipireddy, Ramakrishna

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA