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

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

Mesh adaptation is essential for accurate computational fluid dynamics (CFD) simulations, especially when computational resources are limited.In simulating the advection-diffusion equation using the finite volume method (FVM), mesh adaptation is often necessary to accurately capture sharp gradients and complex flow features.Although the metric-based mesh adaptation used in this study is highly effective with true errors, these are typically unavailable, requiring the use of error estimators.However, adjoint-based error estimators, though effective, are computationally expensive as they require solving the larger adjoint system.To address this, an inexpensive error estimator based on the gradient of the numerical solution and the embedded method is proposed.Specifically, the central difference and upwind difference schemes are utilized for error estimation within this framework.Since errors are predominantly influenced by sharp gradients in the scalar boundary layer, a target functional based on the product of the embedded error and the solution gradient to compute the error estimator is proposed.Consequently, these proposed error estimators generate meshes closely resembling those produced using true error values, effectively resolving the scalar boundary layer.As a result, the L 2 and L norms of the error typically decrease with successive adaptation cycles.Furthermore, the L 1 , L 2 , and L norms of error per element also reduce with successive adaptation cycles.In conclusion, the proposed error estimator facilitates mesh adaptation without needing true errors or computationally expensive adjoint-based error estimators, thus ensuring accurate numerical solutions, especially in scenarios with sharp boundary layers.

This paper employs automatic differentiation (AD) in the compressive sensing-based generalized polynomial chaos (gPC) expansion, which computes a sparse approximation of the Quantity of Interest (QoI) using orthogonal polynomials as basis functions. An earlier approach without AD relies on an iterative procedure to refine the solution by approximating the gradient of the QoI. With AD, the gradient can be accurately evaluated, and a set of basis functions of the gPC expansion associated with new random variables can be efficiently identified. The computational complexity of the algorithm using AD is independent of the number of basis functions, whereas an earlier algorithm had complexity proportional to the square of this number. Our test problems include synthetic problems and a high-dimensional stochastic partial differential equation. With the new basis, the coefficient vector in the gPC expansion is sparser than the original basis. We demonstrate that introducing AD can greatly improve the performance by computing solutions 2 to 10 times faster than an earlier approach. The accuracy of the gPC expansion is also improved; sparse gpC expansions are obtained without iterative refinement, even for high dimensions when an earlier approach fails.

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

Numerical simulations and optimisation methods, such as mesh adaptation, rely on the accurate and inexpensive use of error estimation methods. Adjoint-based error estimation is the most accurate method, and generally the most costly. A strong contributor to this cost is the need to compute a higher resolution adjoint solution. Here, it is proposed to use super-resolution neural networks to super-resolve a fine adjoint solution from a lower-cost coarse adjoint solution: a superAdjoint. The method is compared to reference error estimators on an unsteady Burgers’ equation using the method of manufactured solutions. Two forms of the superAdjoint were implemented, a twice and a four times refining super-resolution neural network. These were used to demonstrate both the computational cost reduction and the potential for the reduction of the storage footprint of the primal problem. The first, referred to as 2×CNN, was able to reconstruct the spatially enriched adjoint solution, thus providing a robust and inexpensive local output error. The second, the 4×CNN, was able to demonstrate the reconstruction ability of super-resolution neural networks for higher upscaling factors. This was leveraged in order to subsample the primal solution in space, thus reducing substantially the storage footprint of the discrete primal solution. Both superAdjoints could achieve the desired level of accuracy when compared to the reconstruction of the refined adjoint-based error estimate. Moreover, superAdjoint was shown to be able to generalise to a new QoI that was untrained for. This gives great confidence in the use of super-resolution neural networks for the reduction of both computational cost and storage requirements of adjoint-based error estimation, and goal-oriented mesh adaptation.

This document summarizes the findings of a strategic planning exercise commissioned by the Weapons Simulation and Computing, Computational Physics (WSC/CP) program at the Lawrence Livermore National Laboratory (LLNL) in FY24. During the year, the committee met with multiple stakeholder communities to gather input, opinions, suggestions and concerns which have been incorporated throughout this document. The key findings from this exercise are summarized: • The development of multiphysics modelling and simulation (mod/sim) codes and software technologies, their deployment on exascale compute platforms, and their broad adoption across the NNSA is a major success of the Advanced Simulation and Computing (ASC) program and the Exascale Computing Project (ECP). Sustained investment in these core technologies is essential. • Today’s state of the art involves running ensembles of O(100K) simulations to perform uncertainty quantification (UQ) and design studies using multiple statistical methods such as Bayesian optimization to understand sensitivities of our models and explore parameterized design spaces. Even with exascale computing, we are practically limited to O(10) parameters in these studies since the number of simulations required to sample the space scales exponentially with the number of design parameters. • The data from these simulation ensembles is increasingly being used to train machine learned (ML) surrogates (or reduced order models, ROMs) which can then be used for optimization or real time design exploration. However, the trained surrogates are still limited in the number of parameters they can represent due to the sampling limitations previously noted. • Augmenting our suite of integrated multiphysics simulation codes, both current and emerging, with the ability to compute gradients (solution derivatives) of arbitrary simulation outputs with respect to (some or all) simulation inputs would be a breakthrough technology, opening the door to a new era of efficient and automated inverse design based on verified and validated mod/sim capabilities. • This capability, which we refer to as differentiable multiphysics codes (DMCs), would revolutionize both UQ and optimization studies by breaking the curse of dimensionality that presently limits our “gradient-free” ensemble based computing approach. A similar breakthrough occurred in the AI/ML community once the ability to compute gradients of arbitrary loss functions using back-propagation became commonplace. Gradient information from the multiphysics codes can also be used to dramatically improve the efficiency and scale of training of ML/ROM surrogates for rapid assessments. • Achieving this in our suite of codes will be a grand challenge, similar to the amount of effort that was required to transition from CPU to GPU computing. It will require buy-in from the entire WSC/CP program and beyond, including all integrated codes, physics and engineering models, third-party library dependencies and performance portability abstractions. It will also require investment in research and development of numerical methods for computing adjoints of coupled physics across multiple adaptively refined moving meshes and of stochastic (Monte Carlo) and mesh free (SPH) methods. • New software and numerical techniques, largely pioneered by the AI/ML community, make this feasible. Chief among these is automatic differentiation (AD), the ability to employ AD at point-wise locations in a physics calculation (instead of traditional black-box approaches) and the ability to perform “back-propagation in time” (or reverse mode AD) for non-linear partial differential equations (PDEs). Fundamentally, the conclusion of this strategic planning exercise is that the time is right to undertake a large scale effort in WSC, centered on the existing integrated codes, to continue the natural evolution of mod/sim in the age of AI/ML. Instead of attempting to replace mod/sim with purely data driven AI/ML models, we believe the key to success is to integrate AI/ML by building on top of the decades of hard-won knowledge and the verified/validated multiphysics modelling capability that is the hallmark of the ASC program.

Error estimation and mesh adaptivity in incompressible viscous flows using a residual power approach

Gradient-based inverse design in photonics has already achieved remarkable results in designing small-footprint, high-performance optical devices. The adjoint variable method, which allows for the efficient computation of gradients, has played a major role in this success. However, gradient-based optimization has not yet been applied to the mode-expansion methods that are the most common approach to studying periodic optical structures like photonic crystals. This is because, in such simulations, the adjoint variable method cannot be defined as explicitly as in standard finite-difference or finite-element time- or frequency-domain methods. Here, we overcome this through the use of automatic differentiation, which is a generalization of the adjoint variable method to arbitrary computational graphs. We implement the plane-wave expansion and the guided-mode expansion methods using an automatic differentiation library, and show that the gradient of any simulation output can be computed efficiently and in parallel with respect to all input parameters. We then use this implementation to optimize the dispersion of a photonic crystal waveguide, and the quality factor of an ultra-small cavity in a lithium niobate slab. This extends photonic inverse design to a whole new class of simulations, and more broadly highlights the importance that automatic differentiation could play in the future for tracking and optimizing complicated physical models.

Advances in Adaptive Computational Methods in Mechanics

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.

Shape sensitivity analysis is a tool that provides quantitative information about the influence of shape parameter changes on the solution of a partial differential equation (PDE). These shape sensitivities are described by a continuous sensitivity equation (CSE). Automatic differentiation (AD) can be used to perform this sensitivity analysis without writing any additional code to solve the sensitivity equation. The approximate solution of the PDE uses a spatial discretization (mesh) that often depends on the shape parameters. Therefore, the straightforward application of AD introduces derivatives of the mesh. There are two drawbacks to this approach. First, extra computational effort (especially memory) is used in these calculations due to mesh sensitivities. Second, this mesh sensitivity information needs to be computed in order to obtain accurate results. In this work, we provide a methodology that avoids mesh sensitivities (and their drawbacks) by defining a modified PDE on a fixed domain (i.e., independent of the shape parameter) such that AD provides the desired approximation of the CSE. Using two examples, we demonstrate significant improvement in the computational effort, both in terms of floating point operations and memory requirements. We explain how these code modifications can be applied to a wide variety of practical problems with minimal changes to the original code. These changes are negligible when compared to the complexity of writing a separate solver for the sensitivity equation.

We review and extend the complex step approach for approximating partial derivatives of smooth vector functions with n independent variables. This method, which leverages complex arithmetic, offers a highly accurate alternative to traditional finite difference techniques, and highly efficient alternative to operator overloading required techniques such as Automatic differentiation. Certain aspects of the multivariable second and higher order derivative approximations are novel. We demonstrate that for 64-bit arithmetic, 15-digit precision is generally obtained for the first partials. 10 to 13 digit precision is routinely obtained for the second partials, and 5 to 8 digit precision is obtained for the third-order partials. The presented formulation does not require operator overloading in modern programming languages such as MATLAB, Python, or C++, and it demands minimal derivation and coding effort. We provide several illustrative examples that demonstrate the utility, accuracy, efficiency, and generality of this approach, including an astrodynamic application.

This work deals with an implementation of automatic differentiation of C++ computer programs in forward mode using operator overloading and expression templates. We report on the efficiency of such implementation and its obvious advantage : the ability to perform sensitivity analysis without touching the source of the computer program by simply adding a library to it. We apply this tool to a flow control problem : minimize the drag of a cylinder, in subsonic unsteady turbulent flow, by controlling the boundary condition of the cylinder.

A number of practical problems in physics can be solved by using accurate higher-order derivatives. Such derivatives can be obtained with automatic differentiation. However, one has to be concerned with the complexity of computing higher-order derivative tensors even for a modest order and number of independents. Initial experiments using univariate Taylor polynomials with interpolation and operator overloading with unrolled loops showed better runtimes than using other automatic differentiation tools. Motivated by these results, we developed the Rapsodia code generator that produces Fortran and C++libraries for the most common intrinsics. Here we explain the algorithmic approach, implementation, and present test results on a select set of applications. Further details on the Rapsodia tool, and an example for user extensions are given in the Appendix.