Programme of the 17th Euro AD Workshop

We revisit an algorithm (called Edge Pushing (EP)) for computing Hessians using Automatic Differentiation (AD) recently proposed by Gower and Mello (2012). Here we give a new, simpler derivation for the EP algorithm based on the notion of live variables from data-flow analysis in compiler theory and redesign the algorithm with close attention to general applicability and performance. We also develop an extension that incorporates preaccumulation to further reduce execution time. Rigorous complexity analyses for the algorithms are provided, and the performance of the algorithms is evaluated using a mesh optimization application and several kinds of synthetic functions as testbeds.

A new model formulation and solution strategy for the design and simulation of chemical processes involving multistream heat exchangers (MHEXs) is presented. The model combines an extension of a classical heat exchanger analysis algorithm with explicit dependence on the heat exchange area in a nonsmooth equation system which can be solved for up to three unknown variables. The method is made tractable by recent advances in the automatic generation of computationally relevant generalized derivatives for piecewise-differentiable functions. A new variant of the vector forward AD mode is applied in conjunction with an extension of the implicit function theorem to calculate these generalized derivatives automatically and efficiently for the MHEX simulation algorithm. A recent LP-Newton method by Facchinei et al. (2014) is then used to solve the nonsmooth equation system with a locally quadratic convergence rate. Our model can efficiently simulate large-scale complex systems, and has the capability to incorporate additional thermodynamic phenomena, such as phase change and equilibrium, which can naturally be described by additional nonsmooth equations.

Piecewise smooth functions (PS) can be approximated locally by piecewise linear models (PL) with a second order error in the distance to the reference point. Those models are locally Lipschitz continuous w.r.t. to the reference point and can be obtained by a slight modification of common AD-Tools as a structured Block matrix called the abs-normal Form. It will be shown under which conditions successive piecewise linearization will lead to a convergent sequence of approximations, whose limit is a root of the underlying PS function.

For the two modes of Automatic Differentiation. The forward mode has been well understood. Arbitrary order of derivatives can be evaluated via propagating a series of Taylor coefficients. In contrast, although first order reverse mode has the nice property that the complexity does not depend on the number of independent variables, how to efficiently evaluate high order derivatives via reverse mode is not well studied. In this work, we introduce the concept of live variable analysis which is used in compiler theory. With the concept, we can give a general loop invariant for reverse mode of AD. By extending the invariant to high orders, we get an algorithm which can evaluate the arbitrary order of derivatives in reverse mode of AD. Further, we give an framework of how to enable the high order reverse mode AD in a MPI program.

One-shot approaches for design optimization augment the solution of the state equation with an adjoint solver yielding approximate derivatives to change the design. The coordination of these iterative processes is well established when only the state equation serves as equality constraint. We propose a modified augmented Lagrangian function for the handling of additional equality constraints. We show that this augmented Lagrangian can be used in gradient-based optimization to solve the original design task.

Julia is a relatively new programming language which has been gaining popularity in many areas of scientific computing. In addition to its main attraction as high-level language which can achieve competitive performance with C and Fortran, Julia is also an ideal platform for developing AD tools which can interact transparently with user code. Our talk will be composed of three sections.

In the first part of the talk, we introduce and interactively demonstrate a number of technical features in Julia, such as type inference, just-in-time compilation, and code introspection which are relevant to implementing AD tools.

In the second part of the talk, we introduce the ForwardDiff.jl package, which implements scalar and vector forward mode of user code. We present recent benchmarks and experiments comparing forward-mode implementations using scalars, stack-allocated vectors, and heap-allocated vectors.

Finally, we present JuMP, a domain-specific modeling language for optimization embedded in Julia. JuMP implements reverse-mode AD on user-provided closed-form expressions to provide sparse Hessians and Jacobians. Efficient computation of the sparse Hessian is enabled by Gebremedhin et al.'s acyclic coloring heuristic. JuMP's AD computations achieve competitive performance with similar commercial tools such as AMPL and GAMS.

The QMUL in-house solver was recently parallelised using OpenMP to exploit parallelism within a compute node in our cluster, and potentially also the XeonPhi cards.

We developed a pre- and postprocessing methodology to maintain the OpenMP parallelism even in the adjoint solver generated with Tapenade. We exploit a symmetry property of our solver algorithm in order to reuse the partitioning from the primal solver in our adjoint solver.

The idea, practical implementation and a case study with run times will be presented.

In this talk we will present some of the ingredients necessary for the application of AD in the context of aerodynamic shape optimization. We discuss some details of the implementation of the AD-based discrete adjoint solver in the open-source framework SU2 and highlight some features of the newly developed CoDiPack library that is an integral part of the solver. Finally, we conclude with some applications for shape optimizations in 3D (turbulent) flows and give an outlook to large-scale applications and optimizations with unsteady flows.

In the problem of optimizing a flow w.r.t. a given set of parameters an alternative to the usage of Reverse Accumulation lies in the utilization of the Discrete Adjoint Formulation. This formulation is obtained by deriving the Jacobian of the residual of the flow equations at a fixed point (the solution of the primal flow equations)[1]. In the context of (steady) CFD simulation both methods have in common that they discard the iteration history and just treat one iteration step with an already converged solution vector. However as opposed to taping the whole iteration step, using the Discrete Adjoint Formulation allows to only consider the calculation of the final residual. I present the implementation and results of this technique in our existing discrete adjoint OpenFOAM Toolkit[2] for a SIMPLE based solver. [1] Discrete Adjoint-Based Approach for Optimization Problems on Three-Dimensional Unstructured Meshes, D. J. Mavriplis, AIAA Journal 2007 [2] A Discrete Adjoint Model for OpenFOAM, M. Towara and U. Naumann, Procedia Computer Science Vol. 18, 2013

Multi-grid methods have been shown to be a very efficient tool for the solution of partial differential equations. Within this talk, we will present the idea of a multi-grid method for non-smooth problems that is based upon techniques from piecewise linear differentiation. In detail, the original problem is approximated by a sequence of piecewise linear models, which can be written in abs-normal form using additional switching variables. In certain cases, one can exploit the structure of the piecewise linearization and formulate the modulus fixed-point iteration for these switching variables in an efficient way. Moreover, one can use the idea of multi-grid methods and find a solution of the modulus fixed-point equation for the switching variables on a coarse discretization, which than serves as an initial guess for the next finer level. Here, the important aspect is the right choice for the prolongation operator to avoid undesirable smoothing effects as it will be shown for some numerical results, which indicate mesh-independent behavior of the resulting method if done in the right way.

Optimization of fluid-flow problems is an active discipline within computational fluid dynamics (CFD), and aims at improving flow devices in terms of e.g. drag reduction, enhancement of mixing, control of laminar-turbulent transition, via the gradient-based optimization of initial/boundary conditions, optimal forcing or optimal shapes.In the context of fully resolved flows, i.e. Direct Numerical Simulations (DNS) of fluid flow phenomena, it is suitable to state the optimization problem as a PDE constrained optimization (here the PDE is given by the Navier-Stokes equations). Using calculus of variations continuous adjoints are computed and subsequently discretized in the same fashion as the PDE itself. Discretizations using a spectral element method are the most accurate choice and lead to minimal accumulations of errors for time dependent problems. We present such a framework as implemented in the code Nek5000 which is meant to allow for a wide class of flow cases and objective functions. Our focus is on nonlinear adjoints which display the added difficulty of frequent checkpointing and thus to slower computation. Using a nearly optimal checkpointing scheme and parallel I/O the compute time of the dual problem can be reduced to 2-5 time the cost of the primal problem. The framework is conceived such that it can handle complex geometries, various boundary conditions and is highly scalable lending itself to large scale simulation. However for validation purposes and efficiency studies we resumed our study to a canonical flow such as the case of a lid driven cavity.

Recently, an interface between ADIC and ADOL-C has been developed at Argonne National Laboratories and the University of Paderborn. This allows applying source transformation techniques to imperative C-style code and, at the same time, operator overloading techniques to object-oriented C++ programs with the aim of improving performance and decreasing memory footprint. On most of these codes ADIC's source transformation cannot be applied at all on its own due to heavy use of object oriented programming and other sophisticated C++ features. However, this potential performance boost needs more empirical verification. As a result of the authors recent stay at Argonne National Laboratory this interface has been applied to the the software packages ISSM, Quantlib and if time allows CFS++. In this talk numerical examples for the potential performance boost are presented.

I want to take the chance to talk about my observations when introducing AD to a python project at Allstate, touch on performance issues in comparison to (py)adolc and the consequences of dynamic typing.

Several deterministic methods for nonconvex optimization require lower bounding information that is typically obtained by minimizing convex underestimators of the objective function on appropriate subdomains. A relaxation scheme by McCormick (1976) is an automatable and efficient method for generating and computing appropriate convex underestimators for composite functions, and is essentially a variant of the forward AD mode. This presentation outlines two new variants of McCormick’s method.

First, a technique for reverse McCormick propagation is presented, to compute convex and concave relaxations for an implicitly-­‐defined set-­‐valued mapping by carrying out certain reverse sweeps of the underlying mapping’s computational graph. This technique can tighten the generated convex underestimators for an objective function of a constrained optimization problem, by exploiting knowledge of the constraints appropriately. The solution rate of a global optimization method using these underestimators can thereby be improved.

Next, modifications are proposed for McCormick’s rules for handling products and compositions with univariate elemental functions. These modified rules eliminate any nonsmoothness introduced through McCormick’s original product and composition rules, and thereby produce convex underestimators that are twice-­‐continuously differentiable, without sacrificing any of the useful properties of McCormick’s method. Gradients may be computed efficiently for these modified underestimators using the standard reverse AD mode, and the underestimators may be minimized using standard smooth optimization techniques.