ylee

The Jacobian and Hessian-free approach for the computational fluid dynamics

Posted on: Mon 26 February 2024

The high-order discrete methods have the potential to increase the solution accuracy in a fixed memory profile. In the field of computational fluid dynamics (CFD), the high-order numerical algorithms need to consider both the spatial and temporal axes, as the solution lies on the spatio-temporal plane. The spatial high-order methods depend on the adjacent data to reconstruct (or interpolate) the solution at the desired point. However, achieving high-order temporal accuracy is more challenging since we don't have any future state information.

One possible way to cope with this challenge is to use the Cauchy-Kowalewski (CK) procedure. With the CK procedure, we can transform the temporal derivatives into the spatial derivatives; thus, we can use them as building blocks for constructing the time Taylor series of the solution. This strategy is known as Lax-Wendroff (LW) type time discretization and is widely studied in the CFD community.

The main benefits of the LW-type method are provided by being a single-stage method. Compared to the multi-stage time integrator, the single-stage methods require fewer data communications between the parallel nodes of the modern high-performance computing architecture. Additionally, the single-stage method is more suitable on the popular grid configuration in the CFD community, Adaptive Mesh Refinement (AMR), where progressively refine the grid resolutions over the simulation runs.

Jacobian and Hessian

Despite the LW-type solvers' advantage, one big implementation hurdle makes them less-portable and less-attractive; they require to compute Jacobian, Hessian, and even more. Although it is possible to compute Jacobian-like tensors with symbolic manipulators like SymPy, Mathematica, and Maple, calculating and implementing them into the code is a cumbersome and error-prone process.

Let us consider the simplest case: 1D Euler equations. This system has three main components: mass conservation, momentum conservation, and energy conservation, so the flux Jacobian is a 3x3 matrix. How about Hessian? Hessian is required to construct third-order temporal accuracy in LW-based schemes, but it is more complicated than the Jacobian matrix; it is rank three tensor. Things are getting worse if we consider multi-dimensional cases. In 3D Euler equations, we need to address 5x5 Jacobian and 5x5x5 Hessian, and they are different in each direction. Of course, for higher than third-order accuracy, higher derivatives of the flux functions are required.

The system-free approach

To meet this end, let me introduce the system-free (SF) approach [1] to bypass the analytical derivations of Jacobian and Hessian; by approximating Jacobian-like tensor-vector contractions. We take the main idea from "Jacobian-Free" Newton-Krylov-type (JFNK) iterative schemes [2], based on the Taylor expansion for the flux vector at a small disturbance from the conservative variables, U \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \mathbf{U}.

F(U+εV)=F(U)+ε∂F∂U⋅V+12ε2∂2F∂U2⋅V⋅V+O(ε3),F(U−εV)=F(U)−ε∂F∂U⋅V+12ε2∂2F∂U2⋅V⋅V+O(ε3). \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \begin{aligned} \mathbf{F} ( \mathbf{U} + \varepsilon \mathbf{V}) &= \mathbf{F}(\mathbf{U}) + \varepsilon \frac{\partial \mathbf{F}}{\partial \mathbf{U}} \cdot \mathbf{V} + \frac{1}{2} \varepsilon^{2} \frac{\partial^{2} \mathbf{F}}{\partial \mathbf{U}^{2}} \cdot \mathbf{V} \cdot \mathbf{V} + \mathcal{O}(\varepsilon^{3}), \\ \mathbf{F} ( \mathbf{U} - \varepsilon \mathbf{V}) &= \mathbf{F}(\mathbf{U}) - \varepsilon \frac{\partial \mathbf{F}}{\partial \mathbf{U}} \cdot \mathbf{V} + \frac{1}{2} \varepsilon^{2} \frac{\partial^{2} \mathbf{F}}{\partial \mathbf{U}^{2}} \cdot \mathbf{V} \cdot \mathbf{V} + \mathcal{O}(\varepsilon^{3}). \end{aligned}

V \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \mathbf{V} is an arbitrary vector with the same number of components as U \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \mathbf{U}, and ε \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \varepsilon is a small scalar perturbation. Note that the dot product (⋅) \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} (\cdot) between Hessian and a vector is to be considered as a tensor contraction. By adding and subtracting the above equations, we would get the following approximations:

∂F∂U⋅V=F(U+εV)−F(U−εV)2ε+O(ε2),∂2F∂U2⋅V⋅V=F(U+εV)−2F(U)+F(U−εV)ε2+O(ε2). \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \begin{aligned} \frac{\partial \mathbf{F}}{\partial \mathbf{U}} \cdot \mathbf{V} &= \frac{\mathbf{F}(\mathbf{U} + \varepsilon \mathbf{V}) - \mathbf{F}(\mathbf{U} - \varepsilon \mathbf{V})}{2\varepsilon} + \mathcal{O}(\varepsilon^{2}), \\ \frac{\partial^{2} \mathbf{F}}{\partial \mathbf{U}^{2}} \cdot \mathbf{V} \cdot \mathbf{V} &= \frac{\mathbf{F}(\mathbf{U} + \varepsilon \mathbf{V}) - 2 \mathbf{F}(\mathbf{U}) + \mathbf{F}(\mathbf{U} - \varepsilon \mathbf{V})}{\varepsilon^{2}} + \mathcal{O}(\varepsilon^{2}). \end{aligned}

If we need to get Hessian-vector-vector contraction with two distinct vectors (V \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \mathbf{V} and W \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \mathbf{W} for example), we can use the simple vector calculus using the Hessian tensor's symmetric property:

∂2F∂U2⋅V⋅W=12(∂2F∂U2⋅(V+W)⋅(V+W)−(∂2F∂U2⋅V⋅V+∂2F∂U2⋅W⋅W))+O(ε2). \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \frac{\partial^{2} \mathbf{F}}{\partial \mathbf{U}^{2}} \cdot \mathbf{V} \cdot \mathbf{W} = \frac{1}{2} \left( \frac{\partial^{2} \mathbf{F}}{\partial \mathbf{U}^{2}} \cdot \left( \mathbf{V} + \mathbf{W} \right) \cdot \left( \mathbf{V} + \mathbf{W} \right) - \left( \frac{\partial^{2} \mathbf{F}}{\partial \mathbf{U}^{2}} \cdot \mathbf{V} \cdot \mathbf{V} + \frac{\partial^{2} \mathbf{F}}{\partial \mathbf{U}^{2}} \cdot \mathbf{W} \cdot \mathbf{W} \right) \right) + \mathcal{O}(\varepsilon^{2}).

Applying the above approximations, we don't need to calculate Jacobian and Hessian directly, but we can approximate them numerically in a more compact form.

Another advantage of the SF procedures is offered by the system independence. Jacobian and Hessian vary with the system of equations, but the SF method is not. Therefore we can utilize the SF method in diverse systems without changing the code too much, from Euler equations to magnetohydrodynamics, for instance.

We have implemented the SF approach to the Picard integration formulation (PIF) method [3], and these works have published in the Journal of Computational Physics. You can find the test results and more detailed informations in our published article.

[1]Lee, Y., & Lee, D. (2021). A single-step third-order temporal discretization with Jacobian-free and Hessian-free formulations for finite difference methods. Journal of Computational Physics, 427, 110063. https://doi.org/10.1016/j.jcp.2020.110063
[2]Knoll, D. A., & Keyes, D. E. (2004). Jacobian-free Newton–Krylov methods: a survey of approaches and applications. Journal of Computational Physics, 193(2), 357-397. https://doi.org/10.1016/j.jcp.2003.08.010
[3]Christlieb, A. J., Guclu, Y., & Seal, D. C. (2015). The Picard integral formulation of weighted essentially nonoscillatory schemes. SIAM Journal on Numerical Analysis, 53(4), 1833-1856. https://doi.org/10.1137/140959936