Algorithmic Differentiation Tool Research Articles

Algorithmic Differentiation for adjoint sensitivity calculation in plasma edge codes

In the framework of fusion energy research, divertor design and model calibration based on plasma edge codes currently rely either on manual iterative tuning of parameters, or alternatively on parameter scans. This, combined with the complex and computationally expensive nature of plasma edge codes, makes these procedures extremely cumbersome. Gradient-based optimization methods can significantly reduce this effort, but require an efficient strategy for sensitivity calculation. Algorithmic Differentiation (AD) offers an efficient and accurate solution, allowing semi-automatic sensitivity computation in complex, continuously developed codes. The adjoint AD mode is especially attractive, as its cost is independent of the number of input parameters. In this paper, adjoint AD sensitivity calculation is deployed for the first time in plasma edge codes, applying the TAPENADE AD tool to SOLPS-ITER. Adjoint AD results are verified to be machine precision accurate compared to tangent AD mode, and up to 10−9 compared to finite differences. Scalings of AD computational efforts prove the advantages of adjoint compared to tangent AD, while memory requirements rapidly increase for adjoint, showing the need for an improved checkpointing strategy. With the new tool, a wealth of information becomes accessible to the modeler. An adjoint sensitivity analysis for a COMPASS density scan identifies the input parameters with largest impact on the solution, and 2D sensitivity maps show their spatial dependence.

Aggregated type handling in CoDiPack

AbstractThe development of algorithmic differentiation (AD) tools focuses mostly on handling floating point types in the target language. Taping optimizations in these tools mostly focus on specific operations like matrix vector products. Aggregated types like std::complex are usually handled by specifying the AD type as a template argument. This approach provides exact results, but prevents the use of expression templates. If AD tools are extended and specialized such that aggregated types can be added to the expression framework, then this will result in reduced memory utilization and improve the timing for applications where aggregated types such as complex number or matrix vector operations are used. Such an integration requires a reformulation of the stored data per expression and a rework of the tape evaluation process. We will demonstrate the overheads on a synthetic benchmark and show the improvement when aggregated types are handled properly by the expression framework of the AD tool.

A novel anisotropic stress‐driven model for bioengineered tissues accounting for remodeling and reorientation based on homeostatic surfaces

AbstractA co‐rotated formulation of the intermediate configuration is derived in a thermodynamically consistent manner. As a result of this formulation, algorithmic differentiation (AD) and the equations of the material model can be combined directly, i.e., the equations can be implemented into the AD tool and the corresponding derivatives can be calculated using AD. This is not possible when the equations are given in terms of the intermediate configuration, since the multiplicative decomposition suffers from an inherent rotational non‐uniqueness. Moreover, a novel stress‐driven kinematic growth model is presented that takes homeostasis and fiber reorientation into account and is based on the co‐rotated formulation. A numerical example reveals the promising potential of both the co‐rotated formulation and the stress‐driven growth model.

BROADCAST: A high-order compressible CFD toolbox for stability and sensitivity using Algorithmic Differentiation

The evolution of any complex dynamical system is described by its state derivative operators. However, the extraction of the exact N-order state derivative operators is often inaccurate and requires approximations. The open-source CFD (Computational Fluid Dynamics) code called BROADCAST discretises the compressible Navier-Stokes equations and then extracts the linearised N-derivative operators through Algorithmic Differentiation (AD) providing a toolbox for laminar flow dynamics. Furthermore, the gradients through adjoint derivation are extracted either by transposition of the linearised operator or through the backward mode of the AD tool. The software includes base-flow computation and linear global stability analysis via eigen-decomposition of the linearised operator or via singular value decomposition of the resolvent operator. Sensitivity tools as well as weakly nonlinear analysis complete the package. The numerical method for the spatial discretisation of the equations consists of a finite-difference high-order shock-capturing scheme applied within a finite volume framework on 2D curvilinear structured grids. The stability and sensitivity tools are demonstrated on two cases: a cylinder flow at low Mach number and a hypersonic boundary layer. Program summaryProgram Title: BROADCASTCPC Library link to program files:https://doi.org/10.17632/hcf893cr2n.1Developer's repository link:https://github.com/OneraHub/BroadcastLicensing provisions: Mozilla Public License 2.0Programming language: Python, FortranNature of problem: The extraction of the exact N-order state derivative operators useful to study complex dynamical systems is often inaccurate and requires approximations from the user.Solution method: This program discretises the compressible Navier-Stokes equations and then extracts the linearised state derivative operators through Algorithmic Differentiation providing a toolbox for laminar flow dynamics.

Hybrid FEM-NN models: Combining artificial neural networks with the finite element method

We present a methodology combining neural networks with physical principle constraints in the form of partial differential equations (PDEs). The approach allows to train neural networks while respecting the PDEs as a strong constraint in the optimisation as apposed to making them part of the loss function. The resulting models are discretised in space by the finite element method (FEM). The method applies to both stationary and transient as well as linear/nonlinear PDEs. We describe implementation of the approach as an extension of the existing FEM framework FEniCS and its algorithmic differentiation tool dolfin-adjoint. Through series of examples we demonstrate capabilities of the approach to recover coefficients and missing PDE operators from observations. Further, the proposed method is compared with alternative methodologies, namely, physics informed neural networks and standard PDE-constrained optimisation. Finally, we demonstrate the method on a complex cardiac cell model problem using deep neural networks.

Adjoint computations by algorithmic differentiation of a parallel solver for time-dependent PDEs

A computational fluid dynamics code is differentiated using algorithmic differentiation (AD) in both tangent and adjoint modes. The two novelties of the present approach are (1) the adjoint code is obtained by letting the AD tool Tapenade invert the complete layer of message passing interface (MPI) communications, and (2) the adjoint code integrates time-dependent, non-linear and dissipative (hence physically irreversible) PDEs with an explicit time integration loop running for ca. 106 time steps. The approach relies on using the Adjoinable MPI library to reverse the non-blocking communication patterns in the original code, and by controlling the memory overhead induced by the time-stepping loop with binomial checkpointing. A description of the necessary code modifications is provided along with the validation of the computed derivatives and a performance comparison of the tangent and adjoint codes.

SICOPOLIS-AD v1: an open-source adjoint modeling framework for ice sheet simulation enabled by the algorithmic differentiation tool OpenAD

Abstract. We present a new capability of the ice sheet model SICOPOLIS that enables flexible adjoint code generation via source transformation using the open-source algorithmic differentiation (AD) tool OpenAD. The adjoint code enables efficient calculation of the sensitivities of a scalar-valued objective function or quantity of interest (QoI) to a range of important, often spatially varying and uncertain model input variables, including initial and boundary conditions, as well as model parameters. Compared to earlier work on the adjoint code generation of SICOPOLIS, our work makes several important advances: (i) it is embedded within the up-to-date trunk of the SICOPOLIS repository – accounting for 1.5 decades of code development and improvements – and is readily available to the wider community; (ii) the AD tool used, OpenAD, is an open-source tool; (iii) the adjoint code developed is applicable to both Greenland and Antarctica, including grounded ice as well as floating ice shelves, with an extended choice of thermodynamical representations. A number of code refactorization steps were required. They are discussed in detail in an Appendix as they hold lessons for the application of AD to legacy codes at large. As an example application, we examine the sensitivity of the total Antarctic Ice Sheet volume to changes in initial ice thickness, austral summer precipitation, and basal and surface temperatures across the ice sheet. Simulations of Antarctica with floating ice shelves show that over 100 years of simulation the sensitivity of total ice sheet volume to the initial ice thickness and precipitation is almost uniformly positive, while the sensitivities to surface and basal temperature are almost uniformly negative. Sensitivity to austral summer precipitation is largest on floating ice shelves from Queen Maud to Queen Mary Land. The largest sensitivity to initial ice thickness is at outlet glaciers around Antarctica. Comparison between total ice sheet volume sensitivities to surface and basal temperature shows that surface temperature sensitivities are higher broadly across the floating ice shelves, while basal temperature sensitivities are highest at the grounding lines of floating ice shelves and outlet glaciers. A uniformly perturbed region of East Antarctica reveals that, among the four control variables tested here, total ice sheet volume is the most sensitive to variations in austral summer precipitation as formulated in SICOPOLIS. Comparison between adjoint- and finite-difference-derived sensitivities shows good agreement, lending confidence that the AD tool is producing correct adjoint code. The new modeling infrastructure is freely available at http://www.sicopolis.net (last access: 2 April 2020) under the development trunk.

Algorithmic differentiation improves the computational efficiency of OpenSim-based trajectory optimization of human movement.

Algorithmic differentiation (AD) is an alternative to finite differences (FD) for evaluating function derivatives. The primary aim of this study was to demonstrate the computational benefits of using AD instead of FD in OpenSim-based trajectory optimization of human movement. The secondary aim was to evaluate computational choices including different AD tools, different linear solvers, and the use of first- or second-order derivatives. First, we enabled the use of AD in OpenSim through a custom source code transformation tool and through the operator overloading tool ADOL-C. Second, we developed an interface between OpenSim and CasADi to solve trajectory optimization problems. Third, we evaluated computational choices through simulations of perturbed balance, two-dimensional predictive simulations of walking, and three-dimensional tracking simulations of walking. We performed all simulations using direct collocation and implicit differential equations. Using AD through our custom tool was between 1.8 ± 0.1 and 17.8 ± 4.9 times faster than using FD, and between 3.6 ± 0.3 and 12.3 ± 1.3 times faster than using AD through ADOL-C. The linear solver efficiency was problem-dependent and no solver was consistently more efficient. Using second-order derivatives was more efficient for balance simulations but less efficient for walking simulations. The walking simulations were physiologically realistic. These results highlight how the use of AD drastically decreases computational time of trajectory optimization problems as compared to more common FD. Overall, combining AD with direct collocation and implicit differential equations decreases the computational burden of trajectory optimization of human movement, which will facilitate their use for biomechanical applications requiring the use of detailed models of the musculoskeletal system.

Discrete adjoint gradient evaluations for linear stress and vibration analysis

This paper presents methods used to perform discrete adjoint gradient evaluations for linear stress and vibration analysis. The methods are implemented within the framework of a discrete adjoint structural solver being developed for multidisciplinary adjoint optimizations of turbomachinery components. The code is differentiated using the algorithmic differentiation tool CoDiPack in tandem with manual treatment of the iterative solvers. Stress analysis leads to a linear system of equations that is typically solved by an iterative solver (e.g. GMRES). To ensure accuracy, the adjoint problem is formulated as a new linear system of equations to be solved. Vibration analysis results in a generalized eigenvalue problem that is also typically solved by an interative solver. The adjoint problem takes out the generalized eigenvalue solve and replaces it by one outer product per eigenfrequency, leading to significantly cheap eigenfrequency gradients for vibration analysis.

Sensitivity analysis of plasma edge code parameters through algorithmic differentiation

Anomalous radial transport coefficients, boundary conditions and reaction rates are among the main sources of uncertainty within plasma edge modeling. In principle, an analysis to determine the sensitivity of code results to changes in uncertain model parameters can be easily implemented through finite differences. However, this incurs in error accumulations and allows scanning only one parameter at a time, requiring a huge computational effort. Algorithmic Differentiation (AD) is a possible alternative to finite differences already applied to several transport codes in different research domains but not yet in plasma edge modeling. AD tools preprocess the source code, identifying elementary operations for which the differential form is well known, and producing a new version of the code that contains the additional derivative information. In this paper, the feasibility of applying AD to plasma edge codes is demonstrated using the TAPENADE tool on the SOLPS-ITER code. As a first preliminary step, the AD tool is applied in the so-called “forward” mode on the B2.5 plasma solver, adopting a fluid neutral approximation to obtain the sensitivities of the calculated quantities of interest on selected code parameters. The proof of principle is carried out by comparing the AD results with those evaluated with finite differences on an ITER H-only case. The sensitivities of the target peak heat load and maximum electron temperature with respect to anomalous radial transport coefficients and core input power are assessed. The comparison with finite differences results in a relative error lower than 10−6. This proves that, in a next step, AD can be exploited for an efficient and accurate sensitivity analysis in the framework of plasma edge simulations.

Expression templates for primal value taping in the reverse mode of algorithmic differentiation

ABSTRACTThe reverse mode of Algorithmic Differentiation (AD) can be implemented in several ways. The major choices are primal value taping vs. Jacobian taping, managed indices vs. unmanaged indices and operator level taping vs. statement level taping. Most of the current AD tools have implemented only one of the eight possible choices, and the data management of the implementation adds another complexity hierarchy. The focus in this paper is the implementation of primal value taping on a statement level. Statement level taping removes the need to create intermediate values on the AD tape which results in reduced memory compared to operator level taping. The implementation will be done for managed and unmanaged indices in the AD tool CoDiPack. Primal value taping with statement level taping has not yet been implemented in any other AD tool, thus we will analyse the properties of the taping approaches and highlight the important details for an efficient implementation. Furthermore, all existing taping approaches in CoDiPack will be compared with the new primal value taping approach. The comparison have been conducted on a simple toy problem and a fully featured computational fluid dynamics solver in the multi-physics suite SU2.

Comparing high-order multivariate AD methods

ABSTRACTTo compute every high-order multivariate derivative value, interpolation methods will be shown to be less accurate than a direct forward multivariate Taylor series method, becoming significant for degrees higher than 10. As order increases, interpolation methods rely on increasingly ill-conditioned matrices where simply rounding exact rational values produced corresponding error in some resulting derivative values. Both interpolation and direct methods use forward AD (algorithmic differentiation); the direct method propagates multivariate series coefficients of the original function, while interpolation methods propagate univariate series of the function in fixed directions and reconstruct the multivariate values. Such interpolation methods, differing in direction choices and reconstruction methods, have been shown to be theoretically more efficient than the direct method for high order. Four alternatives were implemented in MATLAB (interpreted and using random access arrays) on the same laptop. In our implementations, the direct method was competitive and often faster in run time, in addition to maintaining good accuracy. Since AD tools in compiled languages are much faster, more comparison is needed. Direct method efficiency depends on indexing subsets within the large non-rectangular data structure for multivariate series coefficients. We explain key implementation details of our direct method that uses a global reference array.

A benchmark of selected algorithmic differentiation tools on some problems in computer vision and machine learning

Algorithmic differentiation (AD) allows exact computation of derivatives given only an implementation of an objective function. Although many AD tools are available, a proper and efficient implementation of AD methods is not straightforward. The existing tools are often too different to allow for a general test suite. In this paper, we compare 15 ways of computing derivatives including 11 automatic differentiation tools implementing various methods and written in various languages (C++, F#, MATLAB, Julia and Python), 2 symbolic differentiation tools, finite differences and hand-derived computation.We look at three objective functions from computer vision and machine learning. These objectives are for the most part simple, in the sense that no iterative loops are involved, and conditional statements are encapsulated in functions such as abs or logsumexp. However, it is important for the success of AD that such ‘simple’ objective functions are handled efficiently, as so many problems in computer vision and machine learning are of this form.

Adjoint algorithmic differentiation tool support for typical numerical patterns in computational finance

We demonstrate the flexibility and ease of use of C++ algorithmic differentiation (AD) tools based on overloading through application to numerical patterns (kernels) arising in computational finance. While adjoint methods and AD have been known in the finance literature for some time, there are few tools capable of handling and integrating with the C++ codes found in production. Adjoint methods are also known to be very powerful but have potentially infeasible memory requirements. We present several techniques for dealing with this problem and demonstrate them on numerical kernels that occur frequently in finance. We build the discussion around the mature AD tool dco/c++, which is designed to handle arbitrary C++ codes and be highly flexible; however, the sketched concepts can certainly be transferred to other AD solutions including in-house tools. An archive of the source code for the numerical kernels as well as all the AD solutions discussed can be downloaded from an accompanying website. This includes documentation for the code and for dco/c++. Trial licences for dco/c++ are available from Numerical Algorithms Group Ltd.

A usability case study of algorithmic differentiation tools on the ISSM ice sheet model

Algorithmic differentiation (AD) based on operator overloading is often the only feasible approach for applying AD in complex C++ software environments. Challenges pertaining to the introduction of an AD tool based on operator overloading have been studied in the past. However, in order to assess possible performance gains or to verify derivative values, it is advantageous to be able to apply more than one AD tool to a given code. Hence, in this work, we investigate usability issues when exchanging AD tools. Our study is based on the NASA/JPL/UCI Ice Sheet System Model (ISSM) which currently employs the AD tool ADOL-C. We introduce CoDiPack to ISSM, a more recent AD tool offering a similar set of features while promising performance improvements. In addition to the obvious type change for the AD-augmented float type, this transition requires the change to a different adjoint MPI library, adaptation of the MUMPS solver wrapper, and changes to the derivative seeding and extraction routines. We believe that these issues are fairly generic for numerical simulation software, and the issues we report on provide a blueprint for similar undertakings. We also believe that our experiences provide guidance towards the development of AD interfaces that support AD tool interoperability. In addition, we improve upon the memory management of the existing ADOL-C instrumentation, which exhibited considerable runtime problems for higher mesh resolutions. We conduct serial and parallel ISSM model runs on a 2D mass transport benchmark as well as a model of the Pine Island Glacier to verify the derivatives computed by both tools and report on runtime performance and memory usage. In comparison, the CoDiPack AD variant of ISSM runs faster with less memory overhead than the ADOL-C variant and, thus, enables future model runs with an increased number of mesh elements. But the existence of two different AD implementations provides added confidence in the correctness of derivatives, in particular for future AD tool versions.

Source-to-source adjoint Algorithmic Differentiation of an ice sheet model written in C

Algorithmic Differentiation (AD) has become a powerful tool to improve our understanding of the Earth System, because it can generate adjoint code which permits efficient calculation of gradients that are essential to sensitivity studies, inverse problems, parameter estimation and data assimilation. Most source-to-source transformation tools, however, have been designed for FORTRAN and support for C remains limited. Here we use the Adjoinable Land Ice Flow model (ALIF), a C clone of the C++ Ice Sheet System Model (ISSM) and employ source-to-source AD to produce its adjoint code. We present the first running source-to-source adjoint of ALIF, and its application to basal drag inversion under Pine Island Glacier, West Antarctica. ALIF brought several challenges to AD tool development, such as the correct treatment of the context code, which does not compute the differentiable function, but controls this computation through the setup of data structures, including possible aliasing, as well as data-flow reversal in the presence of pointers and dynamic memory, which are ubiquitous in codes such as ISSM and ALIF. We present the strategies we have developed to overcome these challenges.

Sensitivity analysis of the meteorological preprocessor MPP-FMI 3.0 using algorithmic differentiation

Abstract. The meteorological input parameters for urban- and local-scale dispersion models can be evaluated by preprocessing meteorological observations, using a boundary-layer parameterisation model. This study presents a sensitivity analysis of a meteorological preprocessor model (MPP-FMI) that utilises readily available meteorological data as input. The sensitivity of the preprocessor to meteorological input was analysed using algorithmic differentiation (AD). The AD tool used was TAPENADE. The AD method numerically evaluates the partial derivatives of functions that are implemented in a computer program. In this study, we focus on the evaluation of vertical fluxes in the atmosphere and in particular on the sensitivity of the predicted inverse Obukhov length and friction velocity on the model input parameters. The study shows that the estimated inverse Obukhov length and friction velocity are most sensitive to wind speed and second most sensitive to solar irradiation. The dependency on wind speed is most pronounced at low wind speeds. The presented results have implications for improving the meteorological preprocessing models. AD is shown to be an efficient tool for studying the ranges of sensitivities of the predicted parameters on the model input values quantitatively. A wider use of such advanced sensitivity analysis methods could potentially be very useful in analysing and improving the models used in atmospheric sciences.

Reverse-mode algorithmic differentiation of an OpenMP-parallel compressible flow solver

Reverse-mode algorithmic differentiation (AD) is an established method for obtaining adjoint derivatives of computer simulation applications. In computational fluid dynamics (CFD), adjoint derivatives of a cost function output such as drag or lift with respect to design parameters such as surface coordinates or geometry control points are a key ingredient for shape optimization, uncertainty quantification and flow control. The computational cost of CFD applications and their derivatives makes it essential to use high-performance computing hardware efficiently, including multi- and many-core architectures. Nevertheless, OpenMP is not supported in most AD tools, and previously shown methods achieve poor scalability of the derivative code. We present the AD of an OpenMP-parallelized finite volume compressible flow solver for unstructured meshes. Our approach enables us to reuse the parallelization of the original code in the computation of adjoint derivatives. The method works by identifying code segments that can be differentiated in reverse-mode without changing their memory access pattern. The OpenMP parallelization is integrated into the derivative code during the build process in a way that is robust to modifications of the original code and independent of the OpenMP support of the differentiation tool. We show the scalability of our adjoint CFD solver on test cases ranging from thousands to millions of finite volume mesh cells on CPUs with up to 16 threads as well as on an Intel XeonPhi card with 236 threads. We demonstrate that our approach is more practical to implement for production-sized CFD codes and produces more efficient adjoint derivative code than previously shown AD methods.

Algorithmic differentiation for piecewise smooth functions: a case study for robust optimization

This paper presents a minimization method for Lipschitz continuous, piecewise smooth objective functions based on algorithmic differentiation (AD). We assume that all nondifferentiabilities are caused by abs(), min(), and max(). The optimization method generates successively piecewise linearizations in abs-normal form and solves these local subproblems by exploiting the resulting kink structure. Both the generation of the abs-normal form and the exploitation of the kink structure are possible due to extensions of standard AD tools. This work presents corresponding drivers for the AD tool ADOL-C which are embedded in the nonsmooth solver LiPsMin. Finally, minimax problems from robust optimization are considered. Numerical results and a comparison of LiPsMin with other nonsmooth optimization methods are discussed.

Algorithmic Differentiation of Code with Multiple Context-Specific Activities

Algorithmic differentiation (AD) by source-transformation is an established method for computing derivatives of computational algorithms. Static dataflow analysis is commonly used by AD tools to determine the set of active variables, that is, variables that are influenced by the program input in a differentiable way and have a differentiable influence on the program output. In this work, a context-sensitive static analysis combined with procedure cloning is used to generate specialised versions of differentiated procedures for each call site. This enables better detection and elimination of unused computations and memory storage, resulting in performance improvements of the generated code, in both forward- and reverse-mode AD. The implications of this multi-activity AD approach on the static analysis of an AD tool is shown using dataflow equations. The worst-case cost of multi-activity AD on the differentiation process is analysed and practical remedies to avoid running into this worst case are presented. The method was implemented in the AD tool Tapenade, and we present its application to a 3D unstructured compressible flow solver, for which we generate an adjoint solver that performs significantly faster when multi-activity AD is used.