Algorithmic differentiation: Sensitivity analysis and the - PowerPoint PPT Presentation

Algorithmic differentiation: Sensitivity analysis and the computation of adjoints Andrea Walther Institut für Mathematik Universität Paderborn LCCC Workshop on Equation-based Modelling September 19–21, 2012

Outline Introduction Basics of Algorithmic Differentiation (AD) The Forward Mode The Reverse Mode Structure-Exploiting Algorithmic Differentiation Time Structure Exploitation Time and Space Structure Exploitation Conclusions Andrea Walther 1 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Computing Derivatives Sensitivity Optimization Simulation Calculation Andrea Walther 2 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Computing Derivatives Sensitivity Optimization Simulation Calculation input x data modelling computer user program output y theory Andrea Walther 2 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Computing Derivatives Sensitivity Optimization Simulation Calculation input x data modelling computer enhanced user program program output y theory differentiation Andrea Walther 2 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Computing Derivatives Sensitivity Optimization Simulation Calculation input x input x data modelling computer enhanced optimization user user program program algorithm output y output y theory sensitivity ∂ y differentiation ∂ x Andrea Walther 2 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Finite Differences Idea: Taylor-expansion, f : R → R smooth then f ( x + h ) = f ( x ) + hf ′ ( x ) + h 2 f ′′ ( x ) / 2 + h 3 f ′′′ ( x ) / 6 + . . . f ( x + h ) ≈ f ( x ) + hf ′ ( x ) ⇒ Df ( x ) = f ( x + h ) − f ( x ) ⇒ h Andrea Walther 3 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Finite Differences Idea: Taylor-expansion, f : R → R smooth then f ( x + h ) = f ( x ) + hf ′ ( x ) + h 2 f ′′ ( x ) / 2 + h 3 f ′′′ ( x ) / 6 + . . . f ( x + h ) ≈ f ( x ) + hf ′ ( x ) ⇒ Df ( x ) = f ( x + h ) − f ( x ) ⇒ h ◮ simple derivative calculation (only function evaluations!) ◮ inexact derivatives ◮ computation cost often too high F : R n → R OPS ( ∇ F ( x )) ∼ ( n + 1 ) OPS ( F ( x )) ⇒ Andrea Walther 3 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Analytic Differentiation ◮ exact derivatives ◮ f ( x ) = exp ( sin ( x 2 )) ⇒ f ′ ( x ) = exp ( sin ( x 2 )) ∗ cos ( x 2 ) ∗ 2 x Andrea Walther 4 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Analytic Differentiation ◮ exact derivatives ◮ f ( x ) = exp ( sin ( x 2 )) ⇒ f ′ ( x ) = exp ( sin ( x 2 )) ∗ cos ( x 2 ) ∗ 2 x x ′ = f ( x , u ) + IC ◮ min J ( x , u ) such that reduced formulation: J ( x , u ) → � J ( u ) J ′ ( u ) based on symbolic adjoint λ ′ = − f x ( x , u ) ⊤ λ + TC � Andrea Walther 4 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Analytic Differentiation ◮ exact derivatives ◮ f ( x ) = exp ( sin ( x 2 )) ⇒ f ′ ( x ) = exp ( sin ( x 2 )) ∗ cos ( x 2 ) ∗ 2 x x ′ = f ( x , u ) + IC ◮ min J ( x , u ) such that reduced formulation: J ( x , u ) → � J ( u ) J ′ ( u ) based on symbolic adjoint λ ′ = − f x ( x , u ) ⊤ λ + TC � ◮ cost (common subexpression, implementation) ◮ legacy code with large number of lines ⇒ closed form expression not available ◮ consistent derivative information?! Andrea Walther 4 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Gedruckt von Andrea Walther Jan 01, 08 21:46 euler2d.c Seite 29/30 Jan 01, 08 21:46 euler2d.c Seite 30/30 read_input_file(argv[1], &code_control); for (it = 0; it < code_control.nsteps[level]; it++) { double residual; code_control.timestep_type = 0; // calculate timestep size like TAU lift = 0.0; drag = 0.0; // read in CFD mesh read_cfd_mesh(code_control.CFDmesh_name, &gridbase); // calculate actual weight of gradient needed for reconstruction grid[0] = gridbase; if (sum_it+first_step <= code_control.start_2nd_order) weight = 0.0; // remove mesh corner points arizing more than once . . . else if (sum_it+first_step < code_control.full_2nd_order) // e.g. for block structured area and at interface between weight = ( double ) (sum_it+first_step − code_control.start_2nd_order) / // block structured and unstructured area (code_control.full_2nd_order − code_control.start_2nd_o remove_double_points( &gridbase, grid); rder); else // write out mesh in tecplot format weight = 1.0; write_pointdata( name, &(grid[0])); // perform a multigrid cycle on current level // calculate metric of finest grid level mg_cycle(grid+level, &code_control, weight, &residual); /* grid[0].xp[ii][1] += 0.00000001; */ calc_metric(&(grid[0]), &code_control); // if current level is finest level, calculate boundary forces puts(" calc_metric ready "); // (lift and drag) if (level == 0) // create coarse meshes for multigrid, calculate their metric calc_forces(grid, &code_control, &lift, &drag); // and initialze forcing functions to zero for (i = 1; i < code_control.nlevels; i++) // set first l2−residual for normalization, if current cycle is { create_coarse_mesh(&(grid[i−1]), &(grid[i])); // the very first of the computation. init2zero(&(grid[i]), grid[i].force); if ((sum_it + first_step) == 0) } first_residual = (fabs(residual) > 1.0e−10) ? residual: 1.0; puts(" create_coarse_mesh ready "); // print out convergence information to file and standard output // initialize flow field on all grid levels to free stream printf(" IT = %d %20.10e %20.10e %20.10e %4.2f\n ", // quantities sum_it, residual / first_residual, lift, drag, weight); fprintf(conv, " %d %20.10e %20.10e %20.10e\n ", for (i = 0; i < code_control.nlevels; i++) init_field(&(grid[i]), &code_control); sum_it+first_step, residual / first_residual, lift, drag); puts(" init_field ready "); sum_it++; } // if selected read restart file if (code_control.restart == 1) // final time of computation read_restart( " restart ", grid, &code_control, t2 = time(&t2); &first_residual, &first_step); // print out time needed for the time loop // calculate primitive variables for all grid levels and printf (" Zeit : %f\n ", difftime(t2, t1)); // initialize states at the boundary last_step = first_step + code_control.nsteps[0] ; for (i = 0; i < code_control.nlevels; i++) { cons2prim(&(grid[i]), &code_control); fclose(conv); init_bdry_states(&(grid[i])); } // map solution from cell centers to vertices center2point(grid); // open file for writing convergence history conv = fopen(" conv.dat ", " w "); // write out field solution fprintf(conv, " title = convergence\n "); write_eulerdata( " euler.dat ", grid, &code_control); fprintf(conv, " variables = iter, l2res, lift, drag\n "); // write out solution on walls write_surf( " euler−surf.dat ", grid, &code_control); level = 0; printf(" will perform %d steps\n ",code_control.nsteps[level]); // write restart file write_restart( " restart ", grid, &code_control, // starting time of computation first_residual, last_step); t1 = time(&t1); return 0; double lift, drag; } // loop over all multigrid cycles Dienstag Januar 01, 2008 euler2d.c 15/15 Andrea Walther 5 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction The “Hello-World”-Example of AD quay Lighthouse Andrea Walther 6 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction The “Hello-World”-Example of AD quay Lighthouse Andrea Walther 6 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction The “Hello-World”-Example of AD y 2 quay y 2 = γ y 1 ω t y 1 ν Lighthouse Andrea Walther 6 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction The “Hello-World”-Example of AD y 2 quay y 2 = γ y 1 ω t y 1 ν Lighthouse ν tan ( ω t ) y 2 = γ ν tan ( ω t ) y 1 = and γ − tan ( ω t ) γ − tan ( ω t ) Andrea Walther 6 / 24 Algorithmic Differentiation LCCC Workshop 2012

Introduction Evaluation Procedure (Lighthouse) v − 3 = x 1 = ν v − 2 = x 2 = γ v − 1 = x 3 = ω ν tan ( ω t ) = x 4 = t v 0 y 1 = γ − tan ( ω t ) v 1 = v − 1 ∗ v 0 ≡ ϕ 1 ( v − 1 , v 0 ) = tan ( v 1 ) v 2 ≡ ϕ 2 ( v 1 ) v 3 = v − 2 − v 2 ≡ ϕ 3 ( v − 2 , v 2 ) y 2 = γ ν tan ( ω t ) v 4 = v − 3 ∗ v 2 ≡ ϕ 4 ( v − 3 , v 2 ) γ − tan ( ω t ) v 5 = v 4 / v 3 ≡ ϕ 5 ( v 4 , v 3 ) = v 5 ∗ v − 2 ≡ ϕ 6 ( v 5 , v − 2 ) v 6 y 1 = v 5 y 2 = v 6 Andrea Walther 7 / 24 Algorithmic Differentiation LCCC Workshop 2012

The Forward Mode Basics of Algorithmic Differentiation Forward Mode of AD F x y Andrea Walther 8 / 24 Algorithmic Differentiation LCCC Workshop 2012

The Forward Mode Basics of Algorithmic Differentiation Forward Mode of AD F x ( t ) y ( t ) Andrea Walther 8 / 24 Algorithmic Differentiation LCCC Workshop 2012

The Forward Mode Basics of Algorithmic Differentiation Forward Mode of AD ˙ x F x ( t ) y ( t ) Andrea Walther 8 / 24 Algorithmic Differentiation LCCC Workshop 2012

The Forward Mode Basics of Algorithmic Differentiation Forward Mode of AD ˙ x F ˙ y x ( t ) y ( t ) ˙ F Andrea Walther 8 / 24 Algorithmic Differentiation LCCC Workshop 2012