Adjoint-based particle forcing reconstruction and uncertainty quantification

The flow forcing of a single smooth spherical object in uniform flow is well-known and described analytically for low Reynolds and Mach number by the Maxey-Riley equation (Maxey & Riley, 1983). Resolution requirements and the complexity of the interaction between wakes, boundary layers, turbulence and shocks, however, limit the accurate measurement and prediction for flow conditions outside of that range. At higher Reynolds numbers, for example, the flow separates, and can become unstable and the forcing depends on subtle linear and non-linear instabilities that can be sensitive to small changes in the flow conditions. For scenarios with high-Reynolds number effects, for example, simulations of explosively dispersed particle fields originating from a high-speed moving source, basic empirical forcing models have been used (Boiko et al., 1997; Jacobs & Don, 2009; Jacobs et al., 2012). However, it is clear that the accuracy of such models is affected by unsteady forcing (Sen et al., 2018a) and unsteady carrier-phase effects (Sen et al., 2019). In fact, high-speed particle flow has spurred model development over the last decade (Bala; Capecelatro et al., 2014; Jacobs et al., 2012) and led to studies that determine either theoretically, experimentally or through machine learning the particle forcing of single particles affected by moving shocks (Kiselev et al., 2006; Ling et al., 2011; Parmar et al., 2009) and arrays of particle affected by shock diffraction (Sen et al., 2017; Ling et al., 2012; Mehta et al., 2016; Capecelatro et al., 2014; Tryggvason & Dabiri, 2013; Regele et al., 2012, 2014).

In addition to the impact of flow conditions, a range of geometric parameters can affect the particle. For example micro- and meso-scopic, but crucial geometric imperfections and deformation of (condensed) particles’ surfaces impacts on momentum and energy exchange with a carrier flow. A harmonic perturbation of cylindrical shape for example can change the dominant vorticity generation by shear to a baroclinic mechanism (Blanco-Casares & Jacobs, 2022) and yield different forcing. Just these few examples show how many parameters and flow conditions can affect forcing models, and it should not be surprising that an accurate and comprehensive forcing model has eluded the community and realistically is not be feasible for all conditions.

As an alternative to understanding the effect of more flow conditions (free-stream turbulence, free-stream velocity profile, stratification etc.) or particle geometry and the formulation of even more additional and new forcing corrections, we might consider the forcing problem from an inverse perspective and resolve the limitations of modeling and experiment through inference and optimization from observational and simulated data using analytical relations and inputs and outputs of simulation and experiment. Could we, for example, infer the forcing from limited trajectory data? Or might it be possible to assume a random forcing models with a confidence interval and improve this forcing with a limited amount of high-resolution data? To this end an intelligent framework is required that systematically and robustly improves forcing models and associated macro-response protocols with a robust uncertainty quantification using limited, experimental data.

The inference of particle forcing from sparse measurements can be accomplished from an inverse-problem perspective. The dependence of the measurements, either individual particle trajectory data or shape of the cloud distribution, onto the forcing law is highly nonlinear, obtained from solving the physical governing equations. In general, for such nonlinear physics-constrained inverse problems, two different classes of methods are considered in past researches. Data-based methods aim at learning massive amount of data as a distribution, often in reduced-order latent space. Methods such as Generative adversarial networks (GANs) (Hassanaly et al., 2022; Silva et al., 2021a, b; Legler & Janjić, 2022), Auto-encoders (Wang et al., 2022b; Peyron et al., 2021; Amendola et al., 2020) or Long-Short-Term Memory (LSTM) (Wan & Sapsis, 2018) can successfully discover the intrinsic coordinates of the data and take into account the effect of observing by evaluating the probability distribution conditioned on those observations. However, giving the efficiency and practicality of these methods, a thorough physical analysis of the results is often difficult because the foundation is not the governing equation. Model-based methods, on the other hand, initiate from the governing equation and often need to solve the adjoint equations to obtain the gradient of the cost function (Wang et al., 2019c). To maintain consistency with the forward solver and ensure numerical stability during optimization, discrete adjoint operator is often derived and implemented (Wang et al., 2019a). Furthermore, the adjoint fields by solving the adjoint equations gives clear interpretation of the physical meaning (Wang et al., 2022a) as the domain of dependence.

One significant benefit of data-based methods is that uncertainty quantification is straightforward under the probabilistic framework. However, for model-based methods, uncertainty quantification is often much more difficult in terms of the adjoint dynamics. For linear dynamics, Bayesian inference with uncertainty quantification can be easily performed using Gaussian Process Regression (GPR) through the aid of adjoint (Gahungu et al., 2022). For nonlinear dynamics, distribution of the control vector deviates from Gaussian and samples of the probability distribution can be draw from Hamiltonian Monte Carlo (HMC) using the gradient Langevin dyanmics (Yang et al., 2021; Herrmann et al., 2019). A remaining hassle would be the computational cost to solve both the forward and adjoint governing equations repetitively.

To account for parametric, empirical and structural uncertainty fluctuations in forcing models, a probabilistic (stochastic), point-particle and multi-scale perspective as was taken in Jacobs & Udaykumar (2019) and Domínguez-Vázquez et al. (2021). Starting from known empirical models and/or surrogate forcing models from high-resolution simulations (Sen et al., 2018b), that account for Reynolds, Mach number and number density, and other parameters, we can describe the forcing and its dependencies within confidence intervals according to a probability density function. Stochastic dynamic models are formulated that propagate this distribution into the stochastic solutions of trajectories (Domínguez-Vázquez et al., 2021).

In this paper, we reconstruct forcing models using a combination of governing dynamic equations that govern the location and velocity of a particle modeled as a singular point with an unknown forcing and adjoint formulations. We develop a theoretical framework that determines forcing for one-way coupled passive particles, under the assumption that the ambient velocity fields are known.

The remainder of the paper is organized as follows: firstly we focus on point-particle and particle forcing models and derive the optimization framework based utilizing the adjoint dynamics. The forcing reconstruction is then practiced in two test cases, the Arnold–Beltrami–Childress (ABC) flow and a decaying isotropic turbulence flow.

When sparse measurements for the particle trajectory are available, information regarding the forcing is indirect. We consider a special case, where the measurement is only available at the final time $t_{m}$. The cost function in this case can be defined as the difference between the computed and measured particle locations at $t_{m}$, namely,

For general cases where measurements are available for several different $t_{m}$, a simple summation of the individual cost functions will be adopted and a similar procedure is still valid.

To satisfy the governing equations 1 in addition to minimize the cost function 4, we augment the cost function with the Lagrange multiplier $\boldsymbol{x}^{\dagger}_{p},\ \boldsymbol{u}^{\dagger}_{p}$, which results in the function $\mathcal{H}$,

If the forcing parameters $\boldsymbol{\alpha}$ are perturbed by an amount $\delta\boldsymbol{\alpha}$, the change of forcing would result in a deviation of $\boldsymbol{x}_{p}$ and $\boldsymbol{u}_{p}$ from its original trajectory, perturbed by the amount $\delta\boldsymbol{x}_{p}$ and $\delta\boldsymbol{u}_{p}$, respectively. In particular, the perturbation of the ambient velocity $\boldsymbol{u}\big{\rvert}_{\boldsymbol{x}_{p}}$ can be obtained from

In addition, the change of forcing $f(a)$ can be evaluated as,

The perturbation of the cost function can be evaluated as,

Through integration by parts, we can obtain another form,

In case of the single measurement at $t_{m}$, the derivatives of the cost function are,

in backward time axis $\tau=t_{m}-t$. The initial conditions of the adjoint simulations are defined as,

Although the source term of the adjoint equation, $\partial J/\partial\boldsymbol{x}_{p}$ can be taken into account in the initial condition instead, we here retain the most general form so that the adjoint equations apply to a general cost function. Also notice that the adjoint equation involves a singularity around $a=|\boldsymbol{u}-\boldsymbol{u}_{p}|=0$, corresponding to the vanishing forcing when $\boldsymbol{u}-\boldsymbol{u}_{p}=\boldsymbol{0}$ in the forward equation.

By solving the adjoint equations backwards in time, we obtain the final gradients as,

In the current study, the discretization of the adjoint equations are derived from summation by parts using the discretized forward equations. This approach, instead of discretizing the adjoint equations derived from integration by parts, is called the discrete adjoint (Wang et al., 2019b). The discrete adjoint calculates the gradient directions $\partial H/\partial\boldsymbol{\alpha}$ accurately to machine zero and guarantees the convergence of the optimization algorithm.

Hamiltonian Monte Carlo, which is known as a golden approach for sampling from posterior distributions, is an efficient Markov Chain Monte Carlo (MCMC) method based on the Hamiltonian dynamics. Suppose that the posterior distribution of the unknown flow parameters, $\boldsymbol{\alpha}$, follows the distribution,

Given the form of the cost function $J$ being quadratic, we are sampling from the posterior distribution of the flow state given measurements being Gaussian variables. To sample from this given posterior, HMC constructs a Hamiltonian system with a fictitious momentum variable $\boldsymbol{r}$,

and we will sample from the joint distribution

for which we ignore the part for the momentum and obtain the marginal distribution for $\boldsymbol{\alpha}$. The sampling is done by selecting $\boldsymbol{r}$ from its Gaussian distribution and integrate equation (14) using symplectic algorithms such as leapfrog. In addition, Metropolis-Hastings step is incorporated to enhance the efficiency of the sampling (Chib & Greenberg, 1995). Algorithmically, it enables the original adjoint-based gradient descend algorithm (Wang et al., 2019d, b) to extend to a probabilistic framework by adding randomness.

The overall procedure of HMC is shown in algorithm 1. The gradient of the cost function, $\nabla_{\boldsymbol{\alpha}}$ is evaluated from adjoint simulations and using equation (12). The computational cost for solving the adjoint equation is comparable to solving the forward problem and does not depend on the number of unknown parameters, e.g. the size of $\boldsymbol{\alpha}$. In the current study, the parameters for HMC are set considering the efficiency for both the number of equation-solving and convergence behavior of the Monte Carlo sampling.

We first consider the Arnold–Beltrami–Childress (ABC) flow (Jeffreys, 1928), a three-dimensional steady analytical solution of the Euler equation, with the constants of the carrier flow field chosen as $A=\sqrt{3}$, $B=\sqrt{2}$ and $C=1$.

The flow field is shown in figure 1, with particle originated at rest from $\boldsymbol{x}_{p,0}=(\pi,\pi,\pi)^{\top}$ with different Stokes numbers $St=\{1,2,5,10\}$. As a demonstration of the data assimilation and adjoint-based uncertainty quantification, we focus on a simple one-parameter reconstruction, where the forcing function is assumed to be,

where $\alpha$ is the only unknown parameter. Here we fix the Stokes number $St=1$ while bringing the readers to the equivalent effect between changing $\alpha$ and changing $St$.

For $\alpha$ fluctuating in a range $0.95\leq\alpha\leq 1.05$, the trajectory of the particle, $\boldsymbol{x}_{p}$ various slightly, as shown in figure 1$(a)$ in the upper panel. As the true value of $\alpha$ is $\alpha_{T}=1$, the other perturbed values of $\alpha$ lead to mismatched observations at the time $t=t_{m}=8$. As illustrated in equations (9) and (10), the initial position of the adjoint particles are given by the mismatch of the observation, $\boldsymbol{x}^{\dagger}_{p}(t_{m})=\boldsymbol{x}_{m}-\boldsymbol{x}_{p}(t_{m})$, therefore, the trajectory $\boldsymbol{x}^{\dagger}_{p}+\boldsymbol{x}_{p}$ originates all from the same location $\boldsymbol{x}_{m}$ at the measurement time $t_{m}$, and travels backwards in time, forming a collection of adjoint particle tracers shown in figure 1$(a)$ on the bottom panel. We here stress the important properties between the forward and adjoint trajectories: nearby particles separates according to the same Lypunov exponent: a conclusion directly derived form the transposed dynamics of the adjoint (Wang et al., 2022a).

The gradient of the cost function with respect to the parameter $\alpha$ is then given by the forward and adjoint particle velocities, as shown in equation (12). For different tolerance of observation error $\sigma$, the posterior distribution of the parameter $\alpha$ can be estimated through the HMC algorithm described in section 2.2. For this simple example with only one parameter, the PDF is also evaluated either through a direct mapping between $\alpha$ and $\boldsymbol{x}_{p}(t_{m})$ as a validation for the HMC algorithm. Both results are shown in figure 1$(b)$. Results with HMC with $\sigma=1$ shows good agreement with direct mapping.

Results for larger value of $\sigma$ show increasing level of non-linearity: since the observation error follows a Gaussian distribution, the posterior distribution for the parameter would be perfectly Gaussian if the dynamic system is linear, which is the case for small $\sigma$; the distribution would deviate from Gaussian for larger $\sigma$.

As the next step, we consider the more general case where the forcing function is parametrized by a few Fourier modes, as illustrated in equation (3). In the application of ABC flow, where $Re_{p}$ is below 10, a total number of seven Fourier modes is enough to efficiently represent the shape in (2). The particle is initialized at location $x_{p}=y_{p}=z_{p}=2\pi$ at rest. The Reference Reynolds number is $Re_{\infty}=250$, the Stokes number is unity $St=1$, the particle density ratio is $\rho_{p}=500$ and the non-dimensional particle diameter $d_{p}=0.012$. We assume that the initial location of the particle at rest is known, and the observations of the three coordinates of the particle location is available at time $t=8$. We choose the observation of particle position at the non-dimensional time of eight Stokes time units so that the particle trajectory is dominated by an inertial response to the changes in the flow field.

Using the gradient of the cost function, following equation (12), we adopted the Quasi-Newton iterative algorithm to reconstruct the forcing parameters. The initial guess is a constant forcing $f=\alpha_{0}$, i.e., $\alpha_{1}=\alpha_{2}=\ldots=\alpha_{6}=0$. We start from initial guesses with $\alpha_{0}$ ranging from 0 to 4 and arrive at different solution of the forcing that drives the cost function sufficiently small, namely $J<10^{-6}$.

Results of various reconstructed forcing functions are plotted in figure 2$(a)$ in grey colors, darker color represents larger $\alpha_{0}$ in the initial guess and red line marks the true forcing function. The collection of reconstructed forcing demonstrates the ill-posedness of the problem, where infinite number of possible solutions exist. However, the particle trajectories from the reconstructed forcing are plotted in figure 2$(b)$ follow only slightly different paths that coincide at $t=8$ before they further deviate. This also implies that even with large range of function $f$, the form of the forcing can only permit a limited region of arrival for the particle.

To further quantify the uncertainty of the forcing function $f$, the HMC algorithm is adopted. The noise level $\sigma$ is the key to the efficiency of the algorithm. A noise level that is too large would result in unconstrained forcing function with unreasonable range. On the other hand, a noise level that is too small would casue large curvature in the landscape of the cost function and cause difficulty for the sampling algorithm. With the noise level $\sigma=0.01$, the posterior distribution for the coefficients $\alpha_{i}$ of these modes is evaluated using HMC and shown in figure 3.

In order to obtain a reasonable starting point for the HMC algorithm, we start from the true forcing parameters. In real practice when the true forcing is unknown, this could be replaced by performing a few L-BFGS iterations to update the forcing coefficients from the initial guess until we reach the desired level of tolerance. This procedure, in lieu of the traditional burn-in period, is more efficient in term to get near the true solution and serves as a fair starting point for HMC algorithm.

The posterior probability distribution of the forcing is plotted as colored contours in figure 3$(a)$. The mean and standard deviation of the reconstructed forcing are plotted in figure 3$(b)$. The forcing reconstruction is much more accurate in the range $1<Re_{p}<5$ than at other particle Reynolds numbers. Compared with the history of $Re_{p}$ plotted in the background of figure 3$(a)$, it is apparent that the $Re_{p}$ with low uncertainty of the forcing coincide with the high probability of occurrence along the trajectories.

The distribution of the cost function during the HMC sampling is plotted in figure 3$(d)$, where the dashed lines marks the theoretical results of $P(J)\sim\exp(-J/\sigma^{2})$. The mismatch between the HMC and the theoretical prediction is due to the difficulty of sampling near the optimal parameters for high-dimensional problems due to low volume fraction interior of the typical set (Betancourt, 2017).

The algorithm is also tested in three-dimensional homogeneous isotropic turbulence in a cubic domain $\Omega=[0,2\pi]\times[0,2\pi]\times[0,2\pi]$.

The flow fields are obtained with an in-house discontinuous Galerkin compressible DNS solver (Klose et al. (2020) and references therein). For the definition of the initial condition of the flow field we adopted the one used in Jacobs et al. (2005); Domínguez-Vázquez et al. (2022). The particle is initialized at location $(x_{p},\ y_{p},\ z_{p})^{\top}=(\pi,\ \pi,\ 2\pi/3)^{\top}$ with velocity $(u_{p},\ v_{p},\ w_{p})^{\top}=(0,\ 0,\ 1.5)^{\top}$. The Reference Reynolds number is $Re_{\infty}=2,357$, the Stokes number is unity $St=1$, the particle density ratio is $\rho_{p}=250$ and the non-dimensional particle diameter $d_{p}=0.0055$. The observation time here is selected to be $t=6$, where the particle travel for six Stokes time units being the time scales of changes in the flow on the same order than those on the particle phase ($St=1$).

The results of uncertainty quantification using 15 Fourier modes are shown in figure 4. The uncertainty of the forcing function $f$ is the smallest for the $Re_{p}$ with the most occurrence along the particle trajectory. This result is very similar to that in the ABC flow. The underlying reason is that as the particle move faster than the ambient fluid, the momentum effect becomes more and more important. For large $Re_{p}$, inertia plays a vital role and the effect of forcing become nearly negligible, rendering the determination for large $Re_{p}$ much more challenging. The probability of occurrence for different Reynolds number, $P(Re_{p})$ along the trajectories during the HMC sampling are also plotted in figure 3$(a)$.