Non-Intrusive Solution of Stochastic and Parametric Equations Hermann G. Matthies a ıc Giraldi b , Alexander Litvinenko c , Dishi Liu d , and Anthony Nouy b Lo¨ a Institute of Scientific Computing, TU Braunschweig, Brunswick, Germany b ´ Ecole Centrale de Nantes, GeM, Nantes, France c KAUST, Thuwal, Saudi Arabia d Institute of Aerodynamics and Flow Control, DLR, Brunswick, Germany wire@tu-bs.de http://www.wire.tu-bs.de 13 Overview-Barcelona.tex,v 5.3.1.2 2015/01/06 00:23:34 hgm Exp
2 Overview 1. Parametric equations 2. Stochastic model problem 3. “Plain vanilla” Galerkin 4. To be or no to be intrusive 5. Numerical Comparison 6. Galerkin and low-rank tensor approximation 7. Non-intrusive computation 8. Numerical examples TU Braunschweig Institute of Scientific Computing
3 General mathematical setup Consider operator equation, physical system modelled by A : A ( p ; u ) = f ( p ) u ∈ U , f ∈ F , U — space of states, F = U ∗ — dual space of actions / forcings; Operator and rhs depend on parameters p , well posed for all p ∈ P . Iterative solver — convergent for all values of p — iterates for k = 0 , . . . , u ( k +1) ( p ) = S ( p ; u ( k ) ( p ) , R ( p ; u ( k ) ( p )) , with u ( k ) ( p ) → u ∗ ( p ) , where S is one cycle of the solver, and the residuum: R ( u ( k ) ) := R ( p ; u ( k ) ( p )) := f ( p ) − A ( p ; u ( k ) ) . When the residuum vanishes — R ( p ; u ∗ ( p )) = 0 — the mapping S has a fixed point u ∗ ( p ) = S ( p ; u ∗ ( p ) , 0) . TU Braunschweig Institute of Scientific Computing
4 Model stochastic problem Geometry flow = 0 Sources 2 flow out 1.5 Dirichlet b.c. 0 1 0.5 0.5 1 1.5 0 2 Aquifier 2D Model model with stochastic data, p x ∈ G ⊂ R d −∇ · ( κ ( x, p ) ∇ u ( x, p )) = f ( x, ω ) & b.c. , ( κ ( x, p ) ∇ u ( x, p )) · n = g ( x, p ) , x ∈ Γ ⊂ ∂ G , p ∈ P κ stochastic conductivity, f and g stochastic sinks and sources. One p is a realisation of κ, f, g . TU Braunschweig Institute of Scientific Computing
5 Preconditioned residual In the iteration set u ( k +1) = u ( k ) + ∆u ( k ) with ∆u ( k ) := S ( p ; u ( k ) , R ( p ; u ( k ) )) − u ( k ) , and usually P ( ∆u ( k ) ) = R ( p ; u ( k ) ) , so that S ( p ; u ( k ) ) = u ( k ) + P − 1 ( R ( p ; u ( k ) )) . (list of arguments shortened) Here P is some preconditioner, which may depend on p , the iteration counter k , and on the current iterate u ( k ) ; e.g. in Newton’s method P = D u A ( p ; u ( k ) ) . TU Braunschweig Institute of Scientific Computing
6 Iteration Algorithm: Start with some initial guess u (0) k ← 0 while no convergence do Compute ∆u ( k ) ← S ( p ; u ( k ) , R ( p ; u ( k ) )) − u ( k ) u ( k +1) ← u ( k ) + ∆u ( k ) k ← k + 1 end while Uniform contraction: ∀ p, u, v : � S ( p ; u ( p ) , R ( p ; u ( p ))) − S ( p ; v ( p ) , R ( p ; v ( p ))) � U ≤ ̺ ∗ � u ( p ) − v ( p ) � U . with ̺ ∗ < 1 . TU Braunschweig Institute of Scientific Computing
7 Discretisation I Let S ⊂ R P be an appropriate Hilbert space of real functions on P , look for solution in tensor space U := U ⊗ S so that u ( p ) = � ι u ι ς ι ( p ) . Normally, discretise first U by choosing finite-dimensional U N ⊂ U , but results here independent of that. Direct Integration: To compute Quantity of Interest (QoI) � � Q ( u ) = Q ( u ( p ) , p ) µ (d p ) ≈ w z Q ( u ( p z ) , p z ) P z integrand and u ( p z ) have to be computed for all p z : expensive! But decoupled, non-intrusive solves. Want to replace it by proxy / meta model or emulator u ( p ) ≈ u M ( p ) TU Braunschweig Institute of Scientific Computing
8 Discretisation II (Further) discretise U ⊗ S by choosing S M = span { Ψ α ( p ) } ⊂ S to give U ⊗ S M = U M ⊂ U ⊗ S = U . Ansatz u M ( ω ) = � α u α Ψ α ( p ) ∈ U M , Often Ψ α ( p ) = Ψ α ( θ ( p )) , where θ ( p ) = [ . . . , θ ℓ ( p ) , . . . ] are independent. If Ψ α ( θ ) = � ℓ ψ α ℓ ( θ ℓ ) , where α = ( . . . , α ℓ , . . . ) , then S M = � ℓ S M,ℓ allows for higher degree tensors. Simplest computation for Ψ α to be orthogonal (orthonormal), � e.g. in inner product � φ, ϕ � S = P φ ( p ) ϕ ( p ) µ (d p ) How to determine the unknown coefficients u α ? TU Braunschweig Institute of Scientific Computing
9 Solution procedures I Projection of the solution u ( p ) , or of the residuum R ( p ; u ( p )) . Interpolation: Determine the u α by interpolating condition: ! � ∀ p β : u ( p β ) = u M ( p β ) = u α Ψ α ( p β ) . α Simplest when Kronecker- δ property Ψ α ( p β ) = δ α,β satisfied. Solve equation on interpolation points p β — decoupled, non-intrusive solves. Pseudo-spectral projection: Simple as Ψ α are orthonormal. Compute projection inner product (integral) by quadrature, i.e. � � u α = Ψ α ( p ) u ( p ) µ (d p ) ≈ w z Ψ α ( p z ) u ( p z ) , P z solve equation on quadrature points p z — decoupled, non-intrusive solve. TU Braunschweig Institute of Scientific Computing
10 Solution procedures II Mapping u ( · ) �→ u M ( · ) is a projection Π , and to describe a general projection, choose ˆ S M = span { Φ α ( p ) } , projection orthogonal to ˆ S M : ∀ ϕ ∈ ˆ ˆ S ⊥ S M : � ( I − Π ) u, ϕ � = 0 , i.e. M = im( I − Π ) . Approximation properties are determined by S M , stability by ˆ S M . Collocation / Interpolation i.e. solve equation on collocation / interpolation points p β , i.e. Φ β ( p ) = δ ( p − p β ) : ! � R ( p β ; u M ( p β )) = R ( p β ; u α Ψ α ( p β )) = 0 . α With Kronecker- δ : R ( p β ; u β ) = 0 — same as interpolation, decoupled, non-intrusive solve. We worry about the norm � Π � . Norm of collocation projector Π C may grow with M . TU Braunschweig Institute of Scientific Computing
11 Projectors Pseudo-spectral projector Π P is orthogonal, i.e. � Π P � = 1 . This means that ˆ S M = S M , normally Φ α = Ψ α Galerkin : Apply Galerkin weighting. � ∀ β : � Φ β ( p ) , R ( p ; u M ( p )) � = � Φ β ( p ) , R ( p ; u α Ψ α ( p )) � = 0 . α Coupled equations, is it intrusive? When solved in a partitioned way, residua computed by quadrature, it is non-intrusive, needs only residua on qudrature points. To have norm of projector as small as possible (Bubnov-Galerkin), choose orthogonal projection Φ α = Ψ α , TU Braunschweig Institute of Scientific Computing
12 Galerkin on iteration equation α u ( k ) Trick: Project iteration equation. Set u ( k ) ( p ) = � α Ψ α ( p )) and u ( k ) = [ . . . , u ( k ) β , . . . ] : u ( k +1) = u ( k ) + ∆u ( k ) = S ( u ( k ) , R ( u ( k ) )) = ⇒ u ( k +1) = u ( k ) + ∆ M u ( k ) , with ∆ M u ( k ) := [ . . . , � Ψ β , S ( p ; u ( k ) ( p ) , R ( p ; u ( k ) ( p ))) � , . . . ] − u ( k ) Define a mapping S ( u ) : � � S ( u ) := [ . . . , � Ψ β , S ( p ; u α Ψ α ( p ) , R ( p ; u α Ψ α ( p ))) � , . . . ] , α α then ∆ M u ( k ) = S ( u ( k ) ) − u ( k ) and u ( k +1) = u ( k ) + ∆ M u ( k ) = S ( u ( k ) ) . TU Braunschweig Institute of Scientific Computing
13 Convergence Nonlinear block Jacobi algorithm: Start with some initial guess u (0) k ← 0 while no convergence do Compute ∆ M u ( k ) as above u ( k +1) ← u ( k ) + ∆ M u ( k ) k ← k + 1 end while Theorem : The mapping S has the same contraction factor ̺ ∗ . This means that the simple nonlinear block Jacobi algorithm converges as before. TU Braunschweig Institute of Scientific Computing
14 The myth about intrusiveness Folklore: Galerkin methods are intrusive. They can be, but don’t have to. Question: To be or not to be intrusive? Stochastic Galerkin conditions for iteration equation requires S ( u ( k ) ) , approximated by � � � S ( u ( k ) ) ≈ S Z ( u ( k ) ) = p z , u ( k ) ( p z ) , R ( p ; u ( k ) ( p z )) υ z Ψ β ( p z ) S . z to give ∆ M u ( k ) ≈ ∆ Z u ( k ) = S Z ( u ( k ) ) − u ( k ) . This requires the evaluation of preconditioned residuum — one iteration with the solver — at each p z . Theorem still holds with ∆ M u ( k ) replaced by ∆ Z u ( k ) in algorithm. TU Braunschweig Institute of Scientific Computing
15 Numerical example R 4 2 R R 3 R R R R R R 6 1 5 A ( p ; u ) := ( Ku + λ 1 ( p 1 )( u T u ) u ) = λ 2 ( p 2 ) f 0 =: f ( p ) . f 0 := [1 , 0 , 0 , 0 , 0] T . TU Braunschweig Institute of Scientific Computing
16 Numerical example spec Case 1 Case 2 Case 3 Case 4 λ 1 ( p 1 ) p 1 + 2 p 1 + 1 . 1 p 1 + 2 sin(4 p 1 + 2) λ 2 ( p 2 ) p 2 + 25 25 p 2 + 0 . 5 10 p 2 + 30 10 sin( p 2 ) + 30 c.o.v. 2 . 5e − 2 2 . 9e+1 1 . 7e − 1 2 . 2e − 1 1e−02 2nd order polynomial 3rd order polynomial 4th order polynomial 5th order polynomial 1e−04 RMSE 1e−06 1e−08 −10 −8 −6 −4 −2 0 10 10 10 10 10 10 Convergence criteria ( ε tol ) TU Braunschweig Institute of Scientific Computing
Recommend
More recommend
