Automatic differentiation in C++ using expression templates and. application to a flow control problem

Expression templates are a well-known set of techniques for improving the efficiency of operator overloading-based forward mode automatic differentiation schemes in the C\(++\) programming language by translating the differentiation from individual operators to whole expressions. However standard expression template approaches result in a large amount of duplicate computation, particularly for large expression trees, degrading their performance. In this paper we describe several techniques for improving the efficiency of expression templates and their implementation in the automatic differentiation package Sacado (Phipps et al., Advances in automatic differentiation, Lecture notes in computational science and engineering, Springer, Berlin, 2008; Phipps and Gay, Sacado automatic differentiation package. http://trilinos.sandia.gov/packages/sacado/, 2011). We demonstrate their improved efficiency through test functions as well as their application to differentiation of a large-scale fluid dynamics simulation code.

Algorithmic differentiation (AD), also known as automatic differentiation, is a technology for accurate and efficient evaluation of derivatives of a function given as a computer model. The evaluations of such models are essential building blocks in numerous scientific computing and data analysis applications, including optimization, parameter identification, sensitivity analysis, uncertainty quantification, nonlinear equation solving, and integration of differential equations. We provide an introduction to AD and present its basic ideas and techniques, some of its most important results, the implementation paradigms it relies on, the connection it has to other domains including machine learning and parallel computing, and a few of the major open problems in the area. Topics we discuss include: forward mode and reverse mode of AD, higher‐order derivatives, operator overloading and source transformation, sparsity exploitation, checkpointing, cross‐country mode, and differentiating iterative processes.This article is categorized under:Algorithmic Development > Scalable Statistical MethodsTechnologies > Data Preprocessing

In mathematics and computer algebra, automatic differentiation (AD) is a set of techniques to evaluate the derivative of a function specified by a computer program. AD exploits the fact that every computer program, no matter how complicated, executes a sequence of elementary arithmetic operations (addition, subtraction, multiplication, division, etc.), elementary functions (exp, log, sin, cos, etc.) and control flow statements. AD takes source code of a function as input and produces source code of the derived function. By applying the chain rule repeatedly to these operations, derivatives of arbitrary order can be computed automatically, accurately to working precision, and using at most a small constant factor more arithmetic operations than the original program. This paper presents AD techniques available in ROOT, supported by Cling, to produce derivatives of arbitrary C/C++ functions through implementing source code transformation and employing the chain rule of differential calculus in both forward mode and reverse mode. We explain its current integration for gradient computation in TFormula. We demonstrate the correctness and performance improvements in ROOT’s fitting algorithms.

Derivatives play a critical role in computational statistics, examples being Bayesian inference using Hamiltonian Monte Carlo sampling and the training of neural networks. Automatic differentiation (AD) is a powerful tool to automate the calculation of derivatives and is preferable to more traditional methods, especially when differentiating complex algorithms and mathematical functions. The implementation of AD, however, requires some care to insure efficiency. Modern differentiation packages deploy a broad range of computational techniques to improve applicability, run time, and memory management. Among these techniques are operation overloading, region‐based memory, and expression templates. There also exist several mathematical techniques which can yield high performance gains when applied to complex algorithms. For example, semi‐analytical derivatives can reduce by orders of magnitude the runtime required to numerically solve and differentiate an algebraic equation. Open and practical problems include the extension of current packages to provide more specialized routines, and finding optimal methods to perform higher‐order differentiation.This article is categorized under: Algorithmic Development > Scalable Statistical Methods

Sensitivity Derivatives for Advanced CFD Algorithm and Viscous Modeling Parameters via Automatic Differentiation

This paper provides a brief introduction to automatic differentiation and relates it to the tangent linear model and adjoint approaches commonly used in meteorology. After a brief review of the forward and reverse mode of automatic differentiation, the ADIFOR automatic differentiation tool is introduced, and initial results of a sensitivity-enhanced version of the MM5 PSU/NCAR mesoscale weather model are presented. We also present, a novel approach to the computation of gradients that, uses a reverse mode approach at the time loop level and a forward mode approach at every time step. The resulting “pseudoadjoint” shares the characteristic of an adjoint code that the ratio of gradient to function evaluation does not depend on the number of independent variables. In contrast to a true adjoint approach, however, the nonlinearity of the model plays no role in the complexity of the derivative code.Key wordsAutomatic DifferentiationAdjointTangent Linear ModelMM5ADIFORSparsLinCData Assimilation

An efficient implicit high-order flux reconstruction/correction procedure via reconstruction method (FR/CPR) is developed to solve the compressible Navier-Stokes equation on unstructured hexahedral grids. The convergence rate of the proposed solver is greatly improved compared to Lower-upper Symmetry Gauss Seidel (LU-SGS) method. The restart GMRES (restarted Generalized Minimal Residual) is used to tackle the linear system by FR/CPR temporal and spatial discretization, and an automatic differential (AD) by operator overloading is proposed to construct the Jacobian matrix of the discretization system. The effects of automatic differential on the convergence rate of the developed implicit scheme have been investigated. Inviscid and viscous test cases are performed to evaluate the performance of the developed implicit solver.

Automatic differentiation approach for solving one-dimensional flow and heat transfer problems

This article presents our work toward correct and efficient automatic differentiation of OpenMP parallel worksharing loops in forward and reverse mode. Automatic differentiation is a method to obtain gradients of numerical programs, which are crucial in optimization, uncertainty quantification, and machine learning. The computational cost to compute gradients is a common bottleneck in practice. For applications that are parallelized for multicore CPUs or GPUs using OpenMP, one also wishes to compute the gradients in parallel. We propose a framework to reason about the correctness of the generated derivative code, from which we justify our OpenMP extension to the differentiation model. We implement this model in the automatic differentiation tool Tapenade and present test cases that are differentiated following our extended differentiation procedure. Performance of the generated derivative programs in forward and reverse mode is better than sequential, although our reverse mode often scales worse than the input programs.

Automatic differentiation of advanced CFD codes for multidisciplinary design

An introduction to both automatic differentiation and object-oriented programming can enrich a numerical analysis course that typically incorporates numerical differentiation and basic MATLAB computation. Automatic differentiation consists of exact algorithms on floating-point arguments. This implementation overloads standard elementary operators and functions in MATLAB with a derivative rule in addition to the function value; for example, $\sin u$ will also compute $(\cos u)\ast u^{\prime}$, where u and $u^{\prime }$ are numerical values. These methods are mostly one-line programs that operate on a class of value-and-derivative objects, providing a simple example of object-oriented programming in MATLAB using the new (as of release 2008a) class definition structure. The resulting powerful tool computes derivative values and multivariable gradients, and is applied to Newton's method for root-finding in both single and multivariable settings. To compute higher-order derivatives of a single-variable function, another class of series objects keeps Taylor polynomial coefficients up to some order. Overloading multiplication on series objects is a combination (discrete convolution) of coefficients. This idea leads to algorithms for other operations and functions on series objects. A survey of more advanced topics in automatic differentiation includes an introduction to the reverse mode (our implementation is forward mode) and considerations in arbitrary-order multivariable series computation.

For the two-dimensional steady and unsteady turbine cascade flows, the Euler/Navier-Stokes equations with Baldwin-Lomax turbulence model are solved in the Cartesian coordinate system. A locally implicit hybrid algorithm on the mixed type of meshes is employed, where the convective dominated part in the flowfield is studied by TVD scheme to obtain high-resolution results on the triangular elements, and the second- and fourth-order dissipative model is introduced on the O-typed quadrilateral grid in the viscous dominated region to minimize the numerical dissipation. When the steady subsonic and transonic turbulent flows are investigated, the distributions of isentropic Mach number on the blade surface, exit flow angle and loss coefficient are obtained. Comparing the present results with the experimental data, the accuracy and reliability of the current approach is confirmed. By giving a moving wake-type total pressure profile at the inlet plane in the rotor-relative frame of reference, the unsteady transonic inviscid and turbulent flows calculations are performed to understand the interaction of the upstream wake with a moving blade row. The Mach number contours, perturbation component of the unsteady velocity vectors, shear stress and pressure distributions on the blade surface are presented. The physical phenomena, which include periodical flow separation on the suction side, bowing, chopping and distortion of incoming wake, negative jet, convection of the vortices and wake segments, and vortex shedding at the trailing edge, are observed. It is concluded that the unsteady aerodynamic behaviors are strongly dependent on the wake/shock/boundary layer interactions.

We study the correctness of automatic differentiation (AD) in the context of a higher-order, Turing-complete language (PCF with real numbers), both in forward and reverse mode. Our main result is that, under mild hypotheses on the primitive functions included in the language, AD is almost everywhere correct, that is, it computes the derivative or gradient of the program under consideration except for a set of Lebesgue measure zero. Stated otherwise, there are inputs on which AD is incorrect, but the probability of randomly choosing one such input is zero. Our result is in fact more precise, in that the set of failure points admits a more explicit description: for example, in case the primitive functions are just constants, addition and multiplication, the set of points where AD fails is contained in a countable union of zero sets of polynomials.

The unsteady flow in stator–rotor interactions affects the structural integrity, aerodynamic performance of the stages, and blade-surface heat transfer. Numerous viscous and inviscid computer programs are used for the prediction of unsteady flows in two-dimensional and three-dimensional stator–rotor interactions. The relative effects of the various components of flow unsteadiness on heat transfer are under investigation. In this paper it is shown that for subsonic cases, the reduced frequency parameter for boundary-layer calculations is about two orders of magnitude smaller than the reduced frequency parameter for the core flow. This means that for typical stator–rotor interactions, the unsteady flow terms are needed to resolve the location of disturbances in the core flow, but in many cases the instantaneous disturbances can be input in steady-flow boundary-layer computations to evaluate boundary-layer effects in a quasi-steady approximation. This hypothesis is tested by comparing computations with experimental data on a turbine rotor for which there are extensive experimental heat transfer data available in the open literature. An unsteady compressible inviscid two-dimensional computer program is used to predict the propagation of the upstream stator disturbances into the downstream rotor passages. The viscous wake (velocity defect) and potential flow (pressure fluctuation) perturbations from the upstream stator are modeled at the computational rotor–inlet boundary. The effects of these interactions on the unsteady rotor flow result in computed instantaneous velocity and pressure fields. The period of the rotor unsteadiness is one stator pitch. The instantaneous velocity fields on the rotor surfaces are input in a steady-flow differential boundary-layer program, which is used to compute the instantaneous heat transfer rate on the rotor blades. The results of these quasi-steady heat-transfer computations are compared with the results of unsteady heat transfer experiments and with the results of previous unsteady heat transfer computations. The unsteady flow fields explain the unsteady amplitudes and phases of the increases and decreases in instantaneous heat transfer rate. It is concluded that the present method is accurate for quantitative predictions of unsteady heat transfer in subsonic turbine flows.

FDOT: A Fast, memory-efficient and automated approach for Discrete adjoint sensitivity analysis using the Operator overloading Technique