Programme of the 21st Euro AD Workshop

The failure surface of solid materials represents a closed two-dimensional submanifold in three-dimensional principal stress space. For porous ceramics, any stress state which belongs to that surface will lead to a sudden fracture in the bulk or a microstructural collapse. Stress states outside the surface are unattainable. The evaluation of this mechanical characteristic is carried out numerically by the Finite Element Method (FEM) and may be validated by a limited number of mechanical experiments. In the current talk we present a computer code with which we are able to automatically determine elasticity tensors and to automatically construct failure surfaces for complex porous microstructures. The code generates finite element meshes on complex domains and automatically performs a series of finite element simulations with boundary conditions sampling the macroscopic principle stress space. The robust construction of failure surfaces from a finite number of points leads us to a question of automatic differentiation, which we will discuss in some detail. (Joint work with Thomas Hochrainer, Kamen Tushtev, and Kurosch Rezwan.)

In this work we present preliminary results on the use of automatic differentiation (AD) within the scope of radiative transfer problems in atmospheric remote sensing.

The retrieval of atmospheric properties from remote sensing measurements requires to solve an ill-posed problem. Usually the parameters of interest are extracted by minimizing the Tikhonov function using the Gauss-Newton algorithm. Here the radiative transfer forward model is used to compute the top-of-atmosphere radiances for a set of given atmospheric parameters, while the corresponding linearized model is used to build the Jacobian matrix, which consists of partial derivatives of the radiances with respect to atmospheric parameters.

The maintenance of the forward model and its linearized version as two different entities may be a source of issues, because any modification in the radiative transfer forward model must be properly reflected in parallel in its linearized version. The substitution of the linearized forward model with an AD library that overloads mathematical functions is a possible alternative to recycle the forward model code and at the same time to keep the differentiation routines in a separate library.

Starting with a Python-implementation of the radiative transfer model DOME (Discrete Ordinates method with Matrix Exponential) that relies on NumPy, we have developed a library that overloads NumPy functions in a similar way to python-ad, but adding most of the functionality directly to a NumPy ndarray subclass to take advantage of NumPy vectorization. In addition, specific functions required by DOME, e.g. the computation of eigenvalues, have also been overloaded. The Jacobian matrices computed by using AD for test cases have been compared with the results from numerical differentiation in terms of accuracy, time efficiency and memory usage.

Within plasma edge modelling for nuclear fusion application, anomalous radial transport coefficients, boundary conditions and reaction rates are among the main sources of uncertainty. 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) [1] is a possible alternative to finite differences already applied to several transport codes in different research domains and recently applied for the first time in plasma edge modelling [2]. In [2] the feasibility of applying AD to plasma edge codes is demonstrated using the TAPENADE [3] tool on SOLPS-ITER [4], the most recent and advanced version of the SOLPS plasma edge code. The AD tool has been applied in forward mode on the B2.5 plasma solver only to obtain the sensitivities of the calculated quantities of interest on selected code parameters. The proof of principle has been carried out by comparing the AD results with those evaluated with finite differences on a real geometry case, resulting in a relative error lower than 10-6. This has proven that AD can be exploited for an efficient and accurate sensitivity analysis in the framework of plasma edge simulations and, in perspective, for an efficient and reliable tool for optimization studies and parameter estimation, notwithstanding the possibility to make use of the reverse AD.

The new Fortran standards 2003, 2008 and 2018 contain various extensions. For some of them we will present how to handle them in source-to-source AD. Other new features cause issues for AD which can be solved with additional efforts.

Source-Transformation adjoint AD of languages with pointers must restore addresses during the so-called backward sweep. Restoration becomes especially difficult when dynamic memory management has changed the base address of these pointers. Restoration might even prove impossible in languages where pointers are hidden, such as languages with Garbage Collection. We explore the main strategies that may answer this problem, none of which being fully satisfying. We try to exhibit their strengthes and drawbacks, in the hope of helping tools developers pick the one best suited to their differentiation model and targetted audience.

Non-Intrusive Least Squares Shadowing (NILSS) is a promising new approach to compute the derivatives of statistically stationary quantities of interest with respect to input parameters in chaotic dynamical systems. NILSS has recently been shown to produce accurate sensitivities in some chaotic fluid flows where conventional adjoint/tangent equation-based sensitivity computation fails. In this talk, I will describe how the tool OpenAD was introduced into the NILSS algorithm.

Algorithmic Differentiation (AD) can be used to automate the generation of derivatives in arbitrary software projects. This will generate maintainable derivatives, that are always consistent with the computation of the software. If a domain specific language (DSL) is used in a software the state of the art approach is to differentiate the DSL library with the same AD tool. The drawback of this solution is the reduced performance since the compiler is no longer able to optimize the e.g. SIMD operations. We present a new approach which integrates the types and operations of the DSL into the AD tool. It will be an operator overloading tool that is generated from an abstract definition of a DSL. This approach enables the compiler to optimize again e.g. for SIMD operation since all calculations are still performed with the original data types. This will also reduce the required memory for AD since the statements inside the DLS implementation are no longer seen by the AD tool. We will present a first implementation and first performance results.

Computer Aided Design (CAD) tools are extensively used to design industrial components. To incorporate them in a gradient-based optimization loop, one requires the shape sensitivities w.r.t. design parameters of a model to be optimized. Within commercial CAD packages, this information is usually not available and calculated using inaccurate finite differences. To obtain the exact derivatives, algorithmic differentiation can be applied if the sources are available. In this study, the Rolls-Royce in-house airfoil design tool Parablading has been differentiated in the forward mode of AD, using the AD software tools ADOL-C and Tapenade. The differentiated Parablading tool has been coupled with a discrete adjoint CFD solver that is part of the Rolls-Royce in-house HYDRA suite of codes, also produced by algorithmic differentiation. This complete differentiated design chain is applied to gradient-based shape optimization of the TU Berlin TurboLab stator test-case to minimize the total pressure loss and exit angle deviation objectives.

Through this presentation (joint work with Bruno Maugars, Bérenger Berthoul) we want to show the added value of automatic differentiation in aerodynamic flow simulation software. We will also show the code architecture constraints entailed. So, we will briefly present the CFD solver elsA from ONERA, the type of specific problems that we have to deal with (large memory cases generally distributed over computing nodes), the use of AD for our solver, the choice of source-to-source transformation and some results achieved.

TBD

Seismic inversion and imaging are adjoint-based optimization problems that processes up to terabytes of data, regularly exceeding the memory capacity of available computers. Data compression is an effective strategy to reduce this memory requirement by a certain factor, particularly if some loss in accuracy is acceptable. A popular alternative is checkpointing, where data is stored at selected points in time, and values at other times are recomputed as needed from the last stored state. This allows arbitrarily large adjoint computations with limited memory, at the cost of additional recomputations. In this paper we combine compression and checkpointing for the first time to compute a realistic seismic inversion. The combination of checkpointing and compression allows larger adjoint computations compared to using only compression, and reduces the recomputation overhead significantly compared to using only checkpointing.

Automatic differentiation does not always result in optimal performance, for example when applied to BLAS routines, linear solvers, or parallelised stencil computations. In these cases it is beneficial to manually differentiate some core routines, or manually optimise the output of a source-transformation AD tool. However, it is easy to introduce bugs when optimising code. We present a number of test cases where we use formal verification tools to show the equivalence of AD-generated and hand-optimised gradient routines, for arbitrary inputs and activity patterns.

For Monte Carlo approximations of integrals of discontinuous functions, black box algorithmic differentiation does not necessarily produce the correct sensitivities. This problem prominently occurs in computational finance where products such as digital and barrier options have discontinuous payoffs and even vanilla options have payoffs with discontinuous first derivatives. We propose to use the adjoints of a smoothed approximation as a control variate for sensitivities produced by the likelihood ratio method. A smoothing control variate significantly reduces the variance of the sensitivity estimators and allows for more stable parameter calibration.

Expressions involving matrix and tensors are often encountered in the area of machine learning for which derivatives need to be computed.

We present the first algorithmic framework for computing matrix and tensor derivatives that extends seamlessly to higher order derivatives. The framework can be used for symbolic as well as for forward and reverse mode automatic differentiation. Experiments show a speedup between two and three orders of magnitude over state-of-the-art frameworks when evaluating higher order derivatives.

An online interface to our framework can be found at    http://www.MatrixCalculus.org.

We consider a generalization of the trapezoidal rule using piecewise linearization for ODEs with Lipschitz continuous RHS. This talk will cover implementation strategies and results concerning the computational cost.

In this short talk, I will introduce the concepts of a novel Matlab-based Generic Inverse Modeling Toolbox. Instead of parsing an existing or new language and extending it to automatic differentiation it is conceptually based on "bricks" and "braches" in the form of a (nested) Matlab cell array, useful for building the forward model. Each brick carries with it the knowledge of its adjunct operator. The currently existing bricks include a number of complex-valued array based operations important in the field of optics. Some bricks also have "coder" functionality which allows the surrounding framework to compile multiple of these into a seamless single-call Cuda program. A variety of Poisson- and Gaussian-based loss functions for inverse modeling are included. The generic toolbox proved useful for various inverse problems in the field of optics and beyond.

Modelica is a modelling language used heavily in the engineering community and industry to simulate DAE based models. The OpenModelica Compiler is an open source compiler implementation for the Modelica language. This compiler analyses the models and generates C/C++ code for the simulation which can then be executed. In order to utilize AD in these simulations as well as do optimization with such models we devised a novell way of generating the trace structure for evaluation in ADOL-C directly from the OpenModelica Compiler, instead of trying to instrument the Compiler generated C/C++ code with the ADOL-C library to trace the functions at runtime. This talk will give some details on the implemetation and the issues.

In recent years, the performance of computing systems is not longer driven by the clock frequency due to physical limits in power consumption and cooling. Instead, power efficient processing units are integrated into so called multi-core processors allowing for concurrent execution of multiple instructions. Algorithms and software packages need to be parallelized to benefit from the parallel hardware. This also applies for applications using Algorithmic Differentiation. We propose a concept for multithreaded evaluation of Jacobians using Algorithmic Differentiation based on operator overloading and its implementation in ADOL-C. Key aspects of the presented work are parallel data structures for the internal representation of the evaluation procedure (tracing), and the parallelization of the Jacobian driver. Moreover, we will present early results w.r.t. scalability of the approach and the introduced overhead due to parallelization.

Many programs rely on OpenMP as an efficient and reliable way for shared memory parallelization. In this talk concepts and challenges on how to apply algorithmic differentiation in this context will be discussed. Furthermore recent work on how to differentiate such a code with the Open Source AD Tool CoDiPack will be presented.

Vector Extensions like avx2 or AVX512 are a common feature in modern CPU architectures, but are rarely used in Algorithmic Differentiation Tools directly. The overloading tool dco/c++ was modified to achieve better performance in primal, tangent and adjoint versions of numerical kernels. Beside trivial primal and tangent speedups, adjoint versions also profit from a vectorized overloading approach. An adjoint version of an Monte Carlo simulation shows speedups up to 3.1, where the theoretical speedup is 4 using avx2.