Second-order derivatives of the general-purpose finite element package SEPRAN via source transformation

Second-order derivatives are crucial ingredients to a variety of numerical methods. Often, they are difficult to get with numerical differentiation by divided differencing. Automatic differentiation provides an alternative by a program transformation capable of evaluating Jacobians, Hessians, or higher-order derivatives of functions given in the form of computer programs. SEPRAN is a general-purpose finite element package written in Fortran 77 used in various scientific areas ranging from fluid dynamics to structural mechanics to electromagnetism. By transforming SEPRAN twice using the automatic differentiation tool ADIFOR, second-order derivatives are evaluated without truncation error. Numerical experiments are reported in which second-order derivatives of a flow field with respect to an inflow velocity are computed, demonstrating the feasibility of this approach.     

TaylUR, an arbitrary-order diagonal automatic differentiation package for Fortran 95

We present TaylUR, a Fortran 95 module to automatically compute the numerical values of a complex-valued function's derivatives with respect to several variables up to an arbitrary order in each variable, but excluding mixed derivatives. Arithmetic operators and Fortran intrinsics are overloaded to act correctly on objects of a defined type taylor, which encodes a function along with its first few derivatives with respect to the user-defined independent variables. Derivatives of products and composite functions are computed using Leibniz's rule and Faà di Bruno's formula. TaylUR makes heavy use of operator overloading and other Fortran 95 features such as elemental functions.

Program summary

Program title: TaylURCatalogue identifier:ADXR_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADXR_v1_0Program obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions:noneProgramming language:Fortran 95Computer:Any computer with a conforming Fortran 95 compilerOperating system:Any system with a conforming Fortran 95 compilerNo. of lines in distributed program, including test data, etc.:6286No. of bytes in distributed program, including test data, etc:14 994Distribution format:tar.gzNature of problem:Problems that require potentially high orders of derivatives with respect to some variables, such as e.g. expansions of Feynman diagrams in particle masses in perturbative Quantum Field Theory, and which cannot be treated using existing Fortran modules for automatic differentiation [C.W. Straka, ADF95: Tool for automatic differentiation of a FORTRAN code designed for large numbers of independent variables, Comput. Phys. Comm. 168 (2005) 123-139, arXiv:cs.MS/0503014; S. Stamatiadis, R. Prosmiti, S.C. Farantos, auto_deriv: Tool for automatic differentiation of a FORTRAN code, Comput. Phys. Comm. 127 (2000) 343-355].Solution method:Arithmetic operators and Fortran intrinsics are overloaded to act correctly on objects of a defined type taylor, which encodes a function along with its first few derivatives with respect to the user-defined independent variables. Derivatives of products and composite functions are computed using Leibniz's rule and Faà di Bruno's formula.Restrictions:Memory and CPU time constraints may restrict the number of variables and Taylor expansion order that can be achieved. Loss of numerical accuracy due to cancellation may become an issue at very high orders.Unusual features:No mixed higher-order derivatives are computed. The complex conjugation operation assumes all independent variables to be real.Running time:The running time of TaylUR operations depends linearly on the number of variables. Its dependence on the Taylor expansion order varies from linear (for linear operations) through quadratic (for multiplication) to exponential (for elementary function calls).     

Looking for narrow interfaces in automatic differentiation using graph drawing

Automatic differentiation is a powerful technique for evaluating derivatives of functions given in the form of a high-level programming language such as Fortran, C, or C++. This technique is superior, in terms of accuracy, to numerical differentiation because it avoids the truncation error involved in divided difference approximations. In automatic differentiation, the program is treated as a potentially very long composition of elementary functions to which the chain rule of differential calculus is applied over and over again. Because of the associativity of the chain rule, there is room for different strategies computing the same numerical results but whose computational cost may vary significantly. Several strategies exploiting high-level structure of the underlying computer code are known to reduce computational cost as opposed to blindly applying automatic differentiation. An example includes “interface contraction” where one takes advantage of the fact that the number of variables passed between subroutines is small compared with the number of propagated directional derivatives. Unfortunately, these so-called narrow interfaces are not immediately available. The present study investigates the use of the VCG graph drawing tool to recognize narrow interfaces in the computational graph, a certain directed acyclic graph used to represent data dependences of variables in the underlying computer code.     

基于OpenCL的自动微分并行实现及其应用

针对如光束平差这样的大规模优化问题,实现基于OpenCL的并行化自动微分。采用更有效的反向计算模式,实现对多参数函数的导数计算。在OpenCL框架下,主机端完成C/C++形式的函数构建以及基于拓扑排序的计算序列生成,设备端按照计算序列完成函数值以及导数的并行计算。测试结果表明,将实现的自动微分应用于光束平差的雅可比矩阵计算后,相比于采用OpenMP的Ceres Solver,运行速度提高了约3.6倍。     

Term Graphs for Computing Derivatives in Imperative Languages

Automatic differentiation is a technique for the rule-based transformation of a subprogram that computes some mathematical function into a subprogram that computes the derivatives of that function. Automatic differentiation algorithms are typically expressed as operating on a weighted term graph called a linearized computational graph. Constructing this weighted term graph for imperative programming languages such as C/C++ and Fortran introduces several challenges. Alias and definition-use information is needed to construct term graphs for individual statements and then combine them into one graph for a collection of statements. Furthermore, the resulting weighted term graph must be represented in a language-independent fashion to enable the use of AD algorithms in tools for various languages. We describe the construction and representation of weighted term graphs for C/C++ and Fortran, as implemented in the ADIC 2.0 and OpenAD/F tools for automatic differentiation.     

Equivalence of continuum and discrete analytic sensitivity methods for nonlinear differential equations

This paper presents results for the study of equivalence between the total form continuum sensitivity equation (CSE) method and the discrete analytic method of shape design sensitivity analysis. For the discrete analytic method, the sensitivity equations are obtained by taking analytic derivatives of the discretized equilibrium equations with respect to the shape design parameters. For the CSE method, the equilibrium equations are firstly differentiated to form a set of linear continuous sensitivity equations and then discretized for solving the shape sensitivities. The sensitivity equations can be derived by taking the material derivatives (total form) or the partial derivatives (local form) of the equilibrium equations. The total form CSE method is shown for the first time to be equivalent, after finite element discretization, to the discrete analytic method for nonlinear second-order differential equations of a particular form with design dependent loads when they use the same: (1) finite element discretization, (2) numerical integration of element matrices, (3) design velocity fields that are linear with respect to the design variable and (4) shape functions for domain transformation and for design velocity field calculations. The shape sensitivity equations for three-dimensional geometric nonlinear elastic structures and linear potential flow are derived by using both total form CSE and discrete analytic method to show the equivalence of the two methods for these specific examples. The accuracy of shape sensitivity analysis is verified by potential flow around an airfoil and a joined beam with an airfoil under gust load. The results show that analytic sensitivity results are consistent with the complex step results.     

Adaptive finite element methods for shape optimization of linearly elastic structures

Grid adaptive methods combined with means for automatic remeshing are applied to problems in shape optimal design of linearly elastic structures. The quantitative effect of element distortion near the design boundaries is identified in terms of interpolation error associated with the finite element discretization. The grid adaptation is itself formulated as a structural optimization problem, with an objective function that reflects the discretization error. A ‘necessary condition’ from this formulation provides the basis for a computational procedure to predict the modified grid.To avoid the sometimes drastic distortion of the FEM grid that might otherwise occur in conjunction with design change, remeshing must be performed at intermediate stages of the overall solution process. In order to produce results for the optimal shape design without interruption in this process, the computer program combines numerical grid generation and automatic remeshing with the grid adaptation and design change. Results for several shape design problems obtained with the use of grid adaptation are compared to computational results predicted from a fixed grid. Both ‘r-’ and ‘h-adaptation’ are tested.     

CAD-based shape optimisation using adjoint sensitivities

Existing approaches to CAD-based design optimisation using adjoint sensitivities are reviewed and their shortcomings are recalled. An alternative approach is presented which uses the control points of the boundary representation (BRep) as design parameters. The sensitivity of the objective function with respect to the design variables is calculated using automatic differentiation (AD). Results for a 2-D aerofoil are presented.     

面向MPP Fortran 的自动数据分布

自动数据分布是面向大规模并行处理MPP(massively parallel processing)系统程序自动并行化的一项关键技术.数据分布方式直接影响着应用程序在MPP系统上的并行执行性能.本文以MPP Fortran为例,详细探讨了自动数据分布的有关技术,如对准分析、分布方式的产生、静态性能评估和数据重新分布等,并提出了相应的算法.这些算法将在作者研制的面向MPP Fortran 的程序自动并行化工具中实现.     

Sensitivity of optimal shapes of artificial grafts with respect to flow parameters

The difficulties arising in the numerical solution of PDE-constrained shape optimization problems are manifold. Key ingredients are the optimization strategy and the shape deformation method. Furthermore, the robustness of the optimal shape with respect to simulation parameters is of great interest. In this paper, we consider fluid flows described by the incompressible Navier–Stokes equations. Previous studies on artificial bypass grafts indicated the need for specific constitutive models to account for the non-Newtonian nature of blood; in particular, the constitutive model was shown to affect the solution of the shape optimization problem. We employ a shape optimization framework that couples a finite element solver with quasi-Newton-type optimizers and a Bezier spline shape parametrization. To compute derivatives of the optimal shapes with respect to viscosity, we transform the entire optimization framework by combining the automatic differentiation tools Adifor2 and TAPENADE. We demonstrate the impact of the geometry parametrization and of geometric constraints on the optimization outcome. Finally, we employ the transformed framework to compute the sensitivity of the optimal shape of bypass grafts with respect to kinematic viscosity. The resulting sensitivities predict very accurately the influence of viscosity changes on the optimal shape. The proposed methodology provides a powerful tool to further investigate the necessity of intricate constitutive models by taking derivatives with respect to model parameters.     

Concurrent wing and high-lift system aerostructural optimization

A method is presented for concurrent aerostructural optimization of wing planform, airfoil and high lift devices. The optimization is defined to minimize the aircraft fuel consumption for cruise, while satisfying the field performance requirements. A coupled adjoint aerostructural tool, that couples a quasi-three-dimensional aerodynamic analysis method with a finite beam element structural analysis is used for this optimization. The Pressure Difference Rule is implemented in the quasi-three-dimensional analysis and is coupled to the aerostructural analysis tool in order to compute the maximum lift coefficient of an elastic wing. The proposed method is able to compute the maximum wing lift coefficient with reasonable accuracy compared to high-fidelity CFD tools that require much higher computational cost. The coupled aerostructural system is solved using the Newton method. The sensitivities of the outputs of the developed tool with respect to the input variables are computed through combined use of the chain rule of differentiation, automatic differentiation and coupled-adjoint method. The results of a sequential optimization, where the wing shape and high lift device shape are optimized sequentially, is compared to the results of simultaneous wing and high lift device optimization.     

Sensitivity analysis of optimum solutions by updating active constraints: application in aircraft structural design

Often the parameters considered as constants in an optimization problem have some uncertainty and it is interesting to know how the optimum solution is modified when these values are changed. The only way to continue having the optimal solution is to perform a new optimization loop, but this may require a high computational effort if the optimization problem is large. However, there are several procedures to obtain the new optimal design, based on getting the sensitivities of design variables and objective function with respect to a fixed parameter. Most of these methods require obtaining second derivatives which has a significant computational cost. This paper uses the feasible direction-based technique updating the active constraints to obtain the approximate optimum design. This procedure only requires the first derivatives and it is noted that the updating set of active constraints improves the result, making possible a greater fixed parameter variation. This methodology is applied to an example of very common structural optimization problems in technical literature and to a real aircraft structure.     

“To be recorded” analysis in reverse-mode automatic differentiation

The automatic generation of adjoints of mathematical models that are implemented as computer programs is receiving increased attention in the scientific and engineering communities. Reverse-mode automatic differentiation is of particular interest for large-scale optimization problems. It allows the computation of gradients at a small constant multiple of the cost for evaluating the objective function itself, independent of the number of input parameters. Source-to-source transformation tools apply simple differentiation rules to generate adjoint codes based on the adjoint version of every statement. In order to guarantee correctness, certain values that are computed and overwritten in the original program must be made available in the adjoint program. For their determination we introduce a static data-flow analysis called “to be recorded” analysis. Possible overestimation of this set must be kept minimal to get efficient adjoint codes. This efficiency is essential for the applicability of source-to-source transformation tools to real-world applications.     

Domain optimization of a multi-element airfoil using automatic differentiation

Automatic differentiation is used to study the dependence of drag and lift of a four-element airfoil with respect to the position of some flaps in Euler and turbulent Navier–Stokes flow. Aim is the optimization of a cost functional based on drag and lift with respect to the position of the flaps as design parameters. The position is defined by the angle of rotation of the flaps around fixed points. A smooth local grid transformation is presented which retains the neighborhood relations of the original unstructured mesh. The two-dimensional (2D) compressible Navier–Stokes/Euler solver NSC2KE from INRIA and a NAG quasi-Netwon optimization routine serve for the solution of the state equations and the optimization, respectively. The gradient information is supplied by the automatic differentiation tool Adifor by Argonne National Laboratory.     

RIOS: efficient I/O in reverse direction

The reverse mode of automatic differentiation executes the adjoint statements induced by each statement in the original program in the reverse order of the original program flow. This program flow reversal commonly requires storage of information on the control flow of the original program. In addition, intermediate values of variables that are overwritten have to be recorded, as these values may later be needed to compute the partial derivatives of the corresponding statement. The stored information will be accessed in reverse order of being written. This runs contrary to many assumptions made in standard implementations of file systems, operating systems, and input/output (I/O) libraries. A common buffering strategy aimed at speeding up future read requests is to employ read‐ahead. This strategy is useful for accesses in forward direction but is considered to be harmful to the performance of the reverse mode. To increase the performance of the reverse mode, it is also advantageous to interleave computations with the data storage and retrieval operations, which can be achieved using multithreading. To this end, we design and implement a novel software called reverse‐mode I/O stream (RIOS) that is adapted to these particular requirements of the reverse mode. We show the advantages of RIOS in two empirical case studies, an artificially constructed example of a typical I/O pattern in the reverse mode and a real‐world example arising from fluid mechanics, which is studied in Fortran90 and in Matlab where the reverse mode is generated via the automatic differentiation tools Tapenade and ADiMat, respectively. Copyright © 2014 John Wiley & Sons, Ltd.     

基于中弧线-厚度函数的翼型形状解析构造法

复杂翼型几何形状的解析表达对叶片的优化设计有重要的意义,文章研究 了用解析函数构造复杂翼型形状的方法。通过对儒科夫斯基翼型函数的简化,得到用中弧线 -厚度函数表示翼型型线的解析表达式,对式中的相关系数和指数进行重新定义和变换,构 造出包括儒科夫斯基翼型的一般翼型型线的解析表达式;通过进一步分离上、下型线并进行 重新组合的方法可构造出更复杂翼型的形状;再通过增加一个独立的厚度函数项的方法,可 构造出具有光滑尾缘形状的翼型。研究表明,复杂翼型的几何形状可通过有限个参数的解析 函数表达,这些参数不仅具有明确的几何意义,而且使用方便,便于调整翼型的局部形状。 文中给出了用翼型、弦长和扭角函数构造风力机叶片解析函数的应用示例。     

Shape optimization with computational fluid dynamics

In many product design and development applications, computational fluid dynamics CFD has become a useful analytical simulation tool. CFD simulations are quite useful in predicting several response parameters for a given design condition. However, like any analysis tool CFD simulations provide limited insight into the design space and the changes needed to find the optimum design parameters.This paper deals with the shape optimization of fluid flows using CFD and numerical optimization techniques. By integrating a commercial optimization code with a CFD code, a CFD shape optimization tool was developed. To study the effectiveness of the developed tool and its ability to produce results with reasonable CPU time, the shape optimization of an airfoil and S-shaped duct are studied with different numbers of design variables. The developed shape optimization tool along with the optimization and CPU time results are discussed.     

Variational sensitivity analysis and design optimization

Variational method (VM) is employed to derive the co-state equations, boundary (transversality) conditions, and functional sensitivity derivatives. The converged solutions of the state equations together with the steady state solution of the co-state equations are integrated along the domain boundary to uniquely determine the functional sensitivity derivatives with respect to the design function. The application of the variational method to aerodynamic shape optimization problems is demonstrated on internal flow problems at supersonic Mach number range of 1.5. Optimization results for flows with and without shock phenomena are presented. The study shows that while maintaining the accuracy of aerodynamical objective function and constraint within the reasonable range for engineering prediction purposes, variational method provides a substantial gain in computational efficiency, i.e., computer time and memory, when compared with the finite difference computations.     

Coupled adjoint aerostructural wing optimization using quasi-three-dimensional aerodynamic analysis

This paper presents a method for wing aerostructural analysis and optimization, which needs much lower computational costs, while computes the wing drag and structural deformation with a level of accuracy comparable to the higher fidelity CFD and FEM tools. A quasi-three-dimensional aerodynamic solver is developed and connected to a finite beam element model for wing aerostructural optimization. In a quasi-three-dimensional approach an inviscid incompressible vortex lattice method is coupled with a viscous compressible airfoil analysis code for drag prediction of a three dimensional wing. The accuracy of the proposed method for wing drag prediction is validated by comparing its results with the results of a higher fidelity CFD analysis. The wing structural deformation as well as the stress distribution in the wingbox structure is computed using a finite beam element model. The Newton method is used to solve the coupled system. The sensitivities of the outputs, for example the wing drag, with respect to the inputs, for example the wing geometry, is computed by a combined use of the coupled adjoint method, automatic differentiation and the chain rule of differentiation. A gradient based optimization is performed using the proposed tool for minimizing the fuel weight of an A320 class aircraft. The optimization resulted in more than 10 % reduction in the aircraft fuel weight by optimizing the wing planform and airfoils shape as well as the wing internal structure.     

RATSNO—an experiment in software adaptability

Adaptable programs are one of the benefits of structured programming. The adaptability of a program is the degree to which it can be transformed into another program that performs a similar, but slightly different, function. While it is clear that the aims of structured programming are best satisfied by the use of modern programming languages, a great number of programmers must use languages such as Fortran. To alleviate this situation, a number of preprocessors have been introduced that give Fortran a more structured facade. This paper describes an experiment performed to test the adaptability of programs written in RATFOR, one of these preprocessors. Judging from the results, the use of a good preprocessor can significantly increase the adaptability of Fortran programs.