c cientific omputing Automatic Differentiation of Computational Fluid Dynamics Package DROPS M. Petera, J. Willkomm EuroAD Workshop 2007 AACHEN
c cientific Outline omputing • Description of an inverse heat transfer problem • AD of DROPS (developed by A. Reusken and his team in Collaborative Research Centre SFB 540 at RWTH Aachen University) EuroAD Workshop 2007 AACHEN
c cientific Falling film experiment omputing A coupled inverse heat transfer and flow problem for the laminar wavy falling film. : heat flux of the heating foil q h EuroAD Workshop 2007 AACHEN
c cientific Inverse heat flow experiment omputing Γ high-resolution temperature measurements on the boundary T m ( x , t ) 1 q = Γ q : heat flow from to the film c foil 2 EuroAD Workshop 2007 AACHEN
c cientific An inverse heat flow problem omputing Γ Optimization : Estimate parameter vector on the q c 2 Γ using high-resolution measurements on the T m ( x , t ) 1 under the assumption that and are known, T q 0 h − → so that: T ( x , t , q ) T ( x , t ) min Γ c m 1 T ( x , t ; q ) The temperature values come from the solution of c the boundary value problem. Minimize 1 2 = − J ( x , t ; q ) : T ( x , t , q ) T ( x , t ) Γ c c m 2 1 L 2 ∇ with Conjugate Gradient Method , using J ( x , t ; q ) c EuroAD Workshop 2007 AACHEN
c cientific Evaluation of gradients omputing First order derivatives of the function with respect to J ( x , t ; q ) c a parameter vector q c dJ ∇ = ∈ 1 x n J R Q dq c Notation: - Resolution in x, y and z direction of the area Ω n , n , n x y z = + ⋅ + ⋅ + - Number of unknowns n ( n 1 ) ( n 1 ) ( n 1 ) X x y z = + ⋅ + Γ Γ - Number of unknowns on or n ( n 1 ) ( n 1 ) q y z 1 2 = ⋅ + - Dimension of the parameter vector n n ( n 1 ) q Q q t c EuroAD Workshop 2007 AACHEN
c cientific Integration algorithm omputing • Simulation time of is divided into intervals [ 0 t , t ] n f t Ax = • In each time step a linear system of equations b i i of a dimension has to be solved: n solve(A,xi,bi ) X , • The right-hand-side b depends on T ( x , t ) − i 1 i and q ( x , t ) q ( x , t ) − c i 1 c i EuroAD Workshop 2007 AACHEN
c cientific Overview omputing INPUT : Finite-Element-Simulation in DROPS: • A large part of the DROPS code is implemented using C++ templates • Heat flow simulation consists of a loop over time steps • A scalar objective function J with many input parameters ADOL-C ( AD by O ver L oading in C ++) (TU Dresden) • employs adouble data type for recording the calculations on so-called tape • Derivative information calculated by evaluation of the tape For a scalar function J , a reverse mode of AD is well suited due to a constant factor for the additional time complexity. EuroAD Workshop 2007 AACHEN
c cientific “Activating” DROPS (I) omputing Activating variables: (Declaration as the adouble data type) • program parameters (input) • result variables (output) • all intermediate variables (which depend on input and contribute to the calculation of the output variables) std::valarray<double> → std::valarray< adouble > Note: The geometry data of the finite element method (FEM) and the numerical variables of the simulations are clearly separated in DROPS, which allowed for an easy activation of the correct variables. EuroAD Workshop 2007 AACHEN
c cientific “Activating” DROPS (II) omputing Activating generic template functions: • call the function f with parameters of type adouble, • an instance of function f is generated, which operates on adouble arguments template <class T> void function f(DROPS::VectorBaseCL<T> const &v) { T sum=0; for (size_t i = 0, i < v.size(); ++i) sum +=v[i]; } DROPS::VectorBaseCL< adouble > adoubleVector; f(adoubleVector); EuroAD Workshop 2007 AACHEN
c cientific Adapting ADOL-C omputing Problem: • The adouble objects have to be created and destroyed satisfying the LIFO (Last In, First Out) order • DROPS code violates the LIFO condition • Memory leakage Solution: • New memory management algorithm: [Stroustrup, 2000], slightly modified pool allocator EuroAD Workshop 2007 AACHEN
AD of the integration process: c cientific checkpointing omputing Problem: the tape becomes very large by Black-box differentiation Solution: checkpointing (employ revolve library [Griewank, Walther, 2000]) • Decomposition of the simulation in time steps, one n t checkpoint per step • With checkpointing technique the tapes are created n t n , ,... − t 1 and can then be erased O ( n ) • One checkpoint of the size in each time step, X - memory space for the whole tape O ( n ) X EuroAD Workshop 2007 AACHEN
AD of the integration process: c cientific hierarchical approach omputing Problem: potentially incorrect AD differentiation of solve(A,xi,bi) Solution: hierarchical differentiation of solve(A,xi,bi) Advantages : • control over convergence • reduce tape size In reverse mode applied to solve(A,xi,bi): dJ = dJ = for a given vector x calculate b i i dx db i i EuroAD Workshop 2007 AACHEN
c cientific Differentiation of solve(A,xi,bi) omputing Using the theorem about implicit functions: = A ( dx ) db − = i i Ax b 0 − = − i i 1 1 ( dx ) ( db ) A + − = i i ( dA ) x A ( dx ) db 0 i i i dJ dJ = − = A A ( dx ) db 0 i i dx db i i i = T x b A i = dA 0 as the Matrix A is constant. = T T Calculate as the solution of b A b x i i i solve(A, bi^T, xi^T) EuroAD Workshop 2007 AACHEN
c cientific Excluding of solve(A,xi,bi ) omputing In practice: • Stop the taping process before the call of solve(A,xi,bi), and start a new tape afterwards • In DROPS solve(A,xi,bi) is the last call in the time step, so no new tape is needed. Advantages: • The size of the tape does not depend on the number of iteration in the forward CG solver, but is constant for each iteration • solve(A, bi^T, xi^T) is iterated, until the desired relative accuracy is reached. EuroAD Workshop 2007 AACHEN
c cientific Careful activation omputing n : number of recorded operations op n : number of active variables act + The size of the tape: O ( n n ) op act n • should be independent of the number of time steps , n act t since n checkpoints are stored, deactivate the -long n t Q vectors = ⋅ + n n ( n 1 ) Q q t n • deactivate matrix A: decrease the number of act EuroAD Workshop 2007 AACHEN
c cientific Summary omputing • Operator overloading with C++ template based code • Need for eliminating LIFO requirement in ADOL-C • Exploitation of the structure of the program to reduce the run time and storage complexity – checkpointing – hierarchical differentiation – careful activation of variables = • Results: A gradient of a scalar function J with respect to n 140 , 000 Q parameters is computed. An additional computational cost with respect to the simulation requires a factor of 7 for the time complexity and 100MB of additional memory space EuroAD Workshop 2007 AACHEN
c cientific Acknowledgements omputing Our work is funded by SFB 540 Thanks to: Kowarz, Walther TU Dresden Gross, Reusken RWTH Aachen EuroAD Workshop 2007 AACHEN
Recommend
More recommend
