c cientific omputing Automatic Differentiation of Computational - - PowerPoint PPT Presentation

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

Automatic Differentiation of Computational Fluid Dynamics Package DROPS

• M. Petera, J. Willkomm

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

Outline

• 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

• mputing

cientific

Falling film experiment

A coupled inverse heat transfer and flow problem for the laminar wavy falling film.

: heat flux of the heating foil

h

q

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

Inverse heat flow experiment

high-resolution temperature measurements on the boundary : heat flow from to the film

) , ( t x Tm

1

Γ

foil c

q q =

2

Γ

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

An inverse heat flow problem

2

2 1

) , ( ) , , ( 2 1 : ) ; , (

L m c c

t x T q t x T q t x J − =

Γ

The temperature values come from the solution of the boundary value problem.

) ; , (

c

q t x T

Optimization : Estimate parameter vector on the using high-resolution measurements

• n the

under the assumption that and are known, so that:

c

q

h

q

2

Γ ) , ( t x Tm T

1

Γ

min ) , ( ) , , (

1

→ −

Γ

t x T q t x T

m c

Minimize with Conjugate Gradient Method, using

) ; , (

c

q t x J ∇

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

Evaluation of gradients

First order derivatives of the function with respect to a parameter vector

) ; , (

c

q t x J

c

q

Q

n 1 x c

R dq dJ J ∈ = ∇

) 1 ( ) 1 ( ) 1 ( + ⋅ + ⋅ + =

z y x X

n n n n ) 1 ( ) 1 ( + ⋅ + =

z y q

n n n ) 1 (

Q

+ ⋅ =

t q

n n n

z y x

n n n , ,

• Resolution in x, y and z direction of the area Ω
• Number of unknowns
• Number of unknowns on or

1

Γ

2

Γ

• Dimension of the parameter vector

c

q

Notation:

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

Integration algorithm

• Simulation time of

is divided into intervals

• In each time step a linear system of equations
• f a dimension

has to be solved:

• The right-hand-side

depends on

] , [ 0

f

t t

t

n

i i

b Ax = ) , (

1 − i

t x T ) , (

1 − i c

t x q ) , (

i c

t x q

, and

X

n

i

b

solve(A,xi,bi)

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

Overview

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 OverLoading 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

• mputing

cientific

“Activating” DROPS (I)

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

• mputing

cientific

“Activating” DROPS (II)

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

• mputing

cientific

Adapting ADOL-C

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

c

• mputing

cientific

AD of the integration process: checkpointing

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

checkpoint per step

• With checkpointing technique the tapes

are created and can then be erased

• One checkpoint of the size

in each time step,

• memory space for the whole tape

,... ,

1 − t t n

n ) (

X

n O ) (

X

n O

t

n

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

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): for a given vector calculate

AD of the integration process: hierarchical approach

i i

x dx dJ =

i i

b db dJ =

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

Using the theorem about implicit functions: as the Matrix A is constant. Calculate as the solution of solve(A, bi^T, xi^T)

) ( ) ( ) ( = − = − + =

Differentiation of solve(A,xi,bi)

−

i i i i i i i

db dx A db dx A x dA b Ax A db dJ dx dJ A db dx db dx A

i i i i i i

= = =

− − 1 1

) ( ) ( ) (

= dA

T i T i

x b A =

i

b

A b x

T i i =

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

Excluding of solve(A,xi,bi )

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

• mputing

cientific

Careful activation

: number of recorded operations : number of active variables The size of the tape:

• should be independent of the number of time steps ,

since checkpoints are stored, deactivate the -long vectors

• deactivate matrix A: decrease the number of
• p

n

act

n ) (

act

• p

n n O +

t

n

t

n

Q

n

act

n

act

n ) 1 ( + ⋅ =

t q Q

n n n

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

Summary

• 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

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

000 , 140 =

Q

n

EuroAD Workshop 2007

AACHEN

c

• mputing

cientific

Acknowledgements

Our work is funded by SFB 540 Thanks to: Kowarz, Walther TU Dresden Gross, Reusken RWTH Aachen