Multi-physics Simulation Guided Generative Diffusion Models with Applications in Fluid and Heat Dynamics

Denoising Diffusion Implicit Models (DDIMs) use a series of Gaussian noise with increasing variance to corrupt the data, then train DeNNs to gradually reconstruct clean data from corrupted ones. The DDIM consists of multiple steps $t=0,1,2,\cdots,T$ (Ho et al., 2020), each one of which corresponds to a specific level of variance in the Gaussian noise. For clarity, we use $\bm{x}_{t,s}$ to denote the step-$t$ diffusion of the state vector observed at time $s$. It is important to note that $t$ signifies the diffusion step, whereas $s$ indicates the actual time within the dynamic system. Fig. 1 provides a pictorial depiction of this process, illustrating both the diffusion mechanism and the temporal evolution of data.

In DDIM, the forward diffusion process is given as,

where $\bm{z}_{t,s}$ are i.i.d. Gaussian noise vectors and $\beta_{t}$ is a predefined constant that determines the noise variances. Similar to Ho et al. (2020), we introduce notation $\alpha_{t}=1-\beta_{t}$ and $\overline{\alpha}_{t}=\prod_{\tau=1}^{t}\alpha_{\tau}$.

In the continuous limit, the forward diffusion process (1) reduces to a stochastic differential equation (SDE) (Song & Ermon, 2019),

where $\bm{w}_{t}$ is the standard Brownian motion or Weiner process. With a slight abuse of notation, the subscript $t$ denotes a continuous variable in (2) and a discrete variable in (1). In literature, (2) is often called a variance-preserving SDE (Song & Ermon, 2019).

In the task of future physics state prediction, historic observations are often available to practitioners. We use $s_{\texttt{ct}}$ to denote the number of context (already observed) observations, and $\bm{x}_{0,1:s_{\texttt{ct}}}$ as a shorthand notation for the concatenated observed vector $[\bm{x}_{0,1},\bm{x}_{0,2},\cdots,\bm{x}_{0,s_{\texttt{ct}}}]$ (see Time evolution in Fig. 1). Mathematically, the task can be formulated as predicting/sampling the state vector $\bm{x}_{0,s^{*}}$ at a time of interest $s^{*}$ from the predictive distribution $p(\bm{x}_{0,s^{*}}|\bm{x}_{0,1:s_{\texttt{ct}}})$, where $s_{\texttt{ct}}<s^{*}$. The target distribution $p(\bm{x}_{0,s^{*}}|\bm{x}_{0,1:s_{\texttt{ct}}})$ is a conditional distribution of future state vector $\bm{x}_{0,s^{*}}$ given the observed state vectors $\bm{x}_{0,1:s_{\texttt{ct}}}$.

It is important to note that for predictions at multiple future time points, or multiple values of $s^{*}$, practitioners can efficiently apply the DDM framework by simply stacking these state vectors and sampling from the the distribution of stacked vectors. Hence, for simplicity and without any loss of generality, we use $\bm{x}_{0,s^{*}}$ to signify the state vector at any given target time $s^{*}$. Typically in autoregressive tasks, $s^{*}$ is set to $s_{\texttt{ct}}+1$.

To generate high-quality samples, DDMs exploit an existing dataset $\mathcal{D}=\{(\bm{x}_{0,1:s_{\texttt{ct}}}^{(i)},\bm{x}_{0,s^{*}}^{(i)})\}_{i=1}^{N}$ to learn the target conditional distribution, where the superscript $(i)$ is the observation index. Different $i$ denotes different collected evolution trajectories of state vectors. $N$ is the total number of trajectories in the training set.

We briefly describe the training objective of DDIM. By iteratively applying (1), one can show that $\bm{x}_{t,s}$ has the same distribution as $\sqrt{\overline{\alpha}_{t}}\bm{x}_{0,s^{*}}+\sqrt{1-\overline{\alpha}_{t}}\bm{\epsilon}$ where $\bm{\epsilon}$ is a vector whose elements are i.i.d. standard Gaussians. DDIM leverages such fact to train an iterative DeNN $\bm{\epsilon}_{\theta}(\cdot)$ that predicts the noise $\bm{\epsilon}$ from the corrupted sample $\bm{x}_{t,s^{*}}$. More specifically, the training objective of DDIM is,

In (3.1), the diffusion and denoising only happen at the target time $s^{*}$. The objective is to minimize the difference between the noise added to the sample and the noise predicted by the denoising network. It is worth noting that theoretically, the denoising objective (3.1) is related to the score function in statistics: roughly speaking, if the sample size goes to infinity and (3.1) is exactly minimized, the optimal $\bm{\epsilon}_{\theta}^{\star}$ becomes (Song & Ermon, 2019),

where $p\left(\bm{x}_{t,s^{*}}|\bm{x}_{0,1:s_{\texttt{ct}}}\right)$ is the p.d.f. of the random vector $\bm{x}_{t,s^{*}}=\sqrt{\overline{\alpha}_{t}}\bm{x}_{0,s^{*}}+\sqrt{1-\overline{\alpha}_{t}}\bm{\epsilon}$, and the score function is the gradient of the logarithm of the p.d.f.

In the inference stage, DDIM (Song et al., 2020) generates high-quality samples from the approximated score functions. More precisely, DDIM samples $\bm{x}_{T,s^{*}}$ from the standard normal distribution, then applies the denoising network $\bm{\epsilon}_{\theta}$ to iteratively denoise $\bm{x}_{t,s^{*}}$:

for $t$ from $T$ to $1$. The coefficient $\frac{1-\alpha_{t}}{\sqrt{1-\overline{\alpha}_{t}}+\sqrt{\alpha_{t}-\overline{\alpha}_{t}}}$ will converge to $\frac{\beta_{t}}{2\sqrt{1-\overline{\alpha}_{t}}}$ when $\beta_{t}$ is small. We obtain $\bm{x}_{0,s^{*}}$ eventually.

From the perspective of SDEs, if the denoising network is properly trained as in (3.1), the sampling rule (5) is a discretized version of the reverse process ODE (2),

Under the dynamics specified by (6), the random vector $\bm{x}_{0,s^{*}}$ will follow the desired predictive distribution $p(\bm{x}_{0,s^{*}}|\bm{x}_{0,1:s_{\texttt{ct}}})$ if properly initialized. Such connection justifies (5) theoretically (Song et al., 2020).

As discussed, the training objective (3.1) and sampling scheme (5) (or the continuous version (6)) only focus on the statistical patterns in the data and can overlook the physics mechanisms, especially when the size of the training dataset $\mathcal{D}$ is not extremely large. Physics simulations can help alleviate the issue. We assume that simulations can make predictions about the future evolution of the system. The output at time $s$ is $\bm{c}_{1,s}$. We further assume here that the physics simulations are inexpensive, allowing simulation predictions to be obtained for each sample in the training and sampling stages with low latency.

With the physics simulator, we can build an augmented training dataset by combining the training data and inexpensive simulation data at target time $s^{*}$, $\mathcal{D}_{\texttt{aug}}=\{(\bm{x}_{0,s^{*}}^{(i)},\bm{x}_{0,1:s_{\texttt{ct}}}^{(i)},\bm{c}_{1,s^{*}}^{(i)})\}_{i=1}^{N}$. The augmented dataset enables us to train a conditional diffusion network $\bm{\epsilon}_{\theta}(\cdot)$ that takes not only the initial frames $\bm{x}_{0,1:s_{\texttt{ct}}}$ but also the simulation prediction $\bm{c}_{1,s^{*}}$ as its contextual input.

The training objective of the inexpensive physics-conditioned diffusion model thus becomes

Similar to DDIM, we use Monte Carlo approximation to minimize the objective. The pseudocode is presented in Algorithm 2.

Intuitively, the physics context $\bm{c}_{1,s^{*}}$ can bring additional physics knowledge to the model, thus augmenting the authenticity of the prediction. The denoising network $\bm{\epsilon}_{\theta}$ is then trained to absorb such physics knowledge. Such intuition is corroborated by our theoretical analysis in Theorem 4.1 (see Section 4).

Accordingly, the iterative sampling rule becomes,

for $t$ from $T$ to $1$, starting from $\bm{x}_{T,s^{*}}\sim\mathcal{N}(0,\mathbf{I})$.

The sampling rule (3.2) is analogous to that in (5). The only difference is that we augment the input to the DeNN with the predictions from the physics simulations. We also provide a pseudocode of the sampling rule as the choice $1$ of Algorithm 3. Our experiments show that the augmented information can significantly improve the sampling performance.

Clearly, the approach above is simple since simulations are cheap, but what if we have simulations that are expensive? Often, practitioners can obtain physics predictions from more expensive but potentially more accurate physics models. We denote the results from an expensive simulator at time $s^{*}$ as $\bm{c}_{2,s^{*}}$. Then, an augmented dataset of trajectories and simulations is $\{\bm{x}_{0,s^{*}}^{(i)},\bm{x}_{0,1:s_{\texttt{ct}}}^{(i)},\bm{c}_{1,s^{*}}^{(i)},\bm{c}_{2,s^{*}}^{(i)}\}_{i=1}^{N}$. However, due to high computation costs or latency, we assume that ${\bm{c}_{2,s^{*}}^{(i)}}$ may only be available for a subset of $i\in\mathcal{S}_{\text{available}}$. Thus, we cannot directly feed $\bm{c}_{2,s^{*}}$ into the DeNN. This section employs a different strategy to handle the case where $\bm{c}_{2,s^{*}}^{(i)}$ is available. Here, we leverage conditional DDMs to guide the diffusion process by both inexpensive and expensive simulators, namely, $p(\bm{x}_{0,s^{*}}|\bm{c}_{1,s^{*}},\bm{c}_{2,s^{*}})$. Notably, our method exploits the sparsely available $\bm{c}_{2,s^{*}}$, decoupled from the training procedure in Section 3.2, ensuring that the unavailability of $\bm{c}_{2,s^{*}}^{(i)}$ does not affect the training of the denoising network $\bm{\epsilon}_{\theta}$.

We motivate our derivation from the reverse diffusion ODE,

One can see that the conditional score function $\nabla_{\bm{x}_{t,s^{*}}}\log p(\bm{x}_{t,s^{*}}|\bm{x}_{0,1:s_{\texttt{ct}}},\bm{c}_{1,s^{*}},\bm{c}_{2,s^{*}})$ replaces its counterpart in (6) and plays the central role in the reverse diffusion process. Therefore, it suffices to derive an estimate of the conditional score function.

From the Bayes’s rule, we know for any $t\geq 0$,

The first term can be approximated by the denoising network trained by Algorithm 2. We use $\bm{g}(\bm{x}_{0,1:s_{\texttt{ct}}},\bm{c}_{2,s^{*}},\bm{c}_{1,s^{*}},t)$ to denote an estimate of the second term,

In literature, there are multiple ways to construct estimates for $\bm{g}$. When the conditional probability $p(\bm{c}_{2,s^{*}}|\bm{x}_{0,s^{*}},\bm{c}_{1,s^{*}})$ is known, we introduce a conceptually simple and computationally tractable procedure inspired by (Chung et al., 2023). The pseudocode is presented in Algorithm 1.

In Algorithm 1, step 2 uses the Tweedie’s formula (Chung et al., 2023) to produce a point estimate of the sample $\hat{\bm{x}}_{0,s^{*},\theta}$ given a noisy sample $\bm{x}_{t,s^{*}}$. It essentially removes the noise from the noisy sample with the noise estimated by the DeNN $\bm{\epsilon}_{\theta}$. Then, step 3 and 4 use the (scaled) gradient over the clean sample $\hat{\bm{x}}_{0,s^{*},\theta}$ to approximate $\nabla_{\bm{x}_{t,s^{*}}}\log p(\bm{c}_{2,s^{*}}|\bm{x}_{t,s^{*}},\bm{c}_{1,s^{*}})$. Such a procedure is easy to implement and works well in practice. We will defer detailed derivations to supplementary materials.

It is the practitioner’s discretion to choose the instantiation of $p(\bm{c}_{2,s^{*}}|\bm{x}_{0,s^{*}},\bm{c}_{1,s^{*}})$ in Algorithm 1. In principle, the conditional probability model should reflect how the expensive simulation $\bm{c}_{2,s^{*}}$ is related to the observations $\bm{x}_{0,s^{*}}$. We will demonstrate two choices of conditional probability models in the numerical experiments about fumes (5.1) and thermal processes (6.2). In section 3.4, we also present some guidelines for designing the conditional model.

Combining Algorithm 1 with the DDIM discretization of (3.3), a discrete update rule for sampling/inference is,

for $t$ from $T$ to $1$. In (LABEL:eqn:expensivediffusionsample),  (Term 1) corresponds to the $\frac{\beta_{t}}{2}\bm{x}_{t,s^{*}}$ component in (3.3), which is essential in the variance-preserving SDE.  (Term 2) approximates the score function that drives the sample $\bm{x}_{t,s^{*}}$ to high-probability regions predicted by the inexpensive physics simulation $\bm{c}_{1,s^{*}}$. The coefficients are consistent with those in DDIM. Furthermore,  (Term 3) represents the conditioning of the expensive physics simulation $\bm{c}_{2,s^{*}}$ that furnishes additional guidance to $\bm{x}_{t,s^{*}}$.  (Term 3) is the major difference between (LABEL:eqn:expensivediffusionsample) and (3.2). The collective effects of three forces encourage the sample to enter regions where statistical patterns and physics knowledge are congruent.

By initializing $\bm{x}_{T,s^{*}}$ from multiple independent samples from the standard normal distribution and iteratively applying (LABEL:eqn:expensivediffusionsample), one can obtain multiple instances of $\bm{x}_{0,s^{*}}$. The sample variances estimated from these instances provide a straightforward characterization for uncertainty quantification.

To summarize, the pseudocodes of the training and sampling algorithms are presented in Algorithm 2 and Algorithm 3.

Leveraging domain-specific knowledge is crucial for determining the exact form of the conditional probability model. In literature, numerous successful examples of such models exist (Jacobsen et al., 2023; Song et al., 2021b).

In scenarios where different physics simulations are independent, it is reasonable to simplify the conditional probability as $p(\bm{c}_{2,s^{*}}|\bm{x}_{0,s^{*}},\bm{c}_{1,s^{*}})=p(\bm{c}_{2,s^{*}}|\bm{x}_{0,s^{*}})$. Our studies demonstrate that this probability can be effectively represented by energy-based models: $p(\bm{c}_{2,s^{*}}|\bm{x}_{0,s^{*}})\propto\exp(-\gamma E(\bm{c}_{2,s^{*}},\bm{x}_{0,s^{*}}))$, where $E$ denotes a differentiable energy function and the partition function is neglected as after taking the logarithm, the partition function contributes only a constant term, which becomes zero if we take gradient on $\bm{x}_{0,s^{*}}$. This energy function attains lower values for consistent pairs of simulation $\bm{c}_{2,s^{*}}$ and sample $\bm{x}_{0,s^{*}}$ and higher values for inconsistent pairs. For instance, if $\bm{c}_{2,s^{*}}$ is a coarse-grained prediction for $\bm{x}_{0,s^{*}}$, one could define $E$ as $E(\bm{c}_{2,s^{*}},\bm{x}_{0,s^{*}})=\left\lVert\bm{c}_{2,s^{*}}-\bm{x}_{0,s^{*}}\right\rVert^{2}$. $\gamma$ serves as a temperature parameter that modulates the strength of the conditional probability. This model framework promotes consistency between the sample and the corresponding expensive physical simulation. We will explore two applications of the energy-based approach in Sections 5 and 6.

In the numerical study, we analyze the movement of 2D fumes driven by buoyancy and gravity, with the goal of predicting buoyancy fields. We use Boussinesq approximation (Spiegel & Veronis, 1960) to analyze the evolution of the fume system. The ground truth data are generated by running multiple simulations of the fumes with existing high-performance CFD simulators (Mohanan et al., 2019b, a) on fine-grained $128\times 128$ grids. More specifically, we randomly initialize the buoyancy and vorticity at time $s=0$ and run the simulator from $s=0$ to $s=10$. We use the buoyancy field at $s=0$ time steps as the context vectors, and predict the target state at $s^{*}=10$. Results of $N=2384$ simulations from different random initializations are accumulated and then randomly separated into 90% training set and 10% test set. We train the DeNN $\bm{\epsilon}_{\theta}$ on the training set. On the test set, we try to predict $\bm{x}_{0,s^{*}}$ with the information $\bm{x}_{0,1:s_{\texttt{ct}}}$, $\bm{c}_{1,s^{*}}$, and possible $\bm{c}_{2,s^{*}}$. The buoyancy at $s^{*}=10$ for four samples from the test set is plotted in the last column of Fig. 2.

To apply MPDM, inexpensive physics predictions $\bm{c}_{1,s^{*}}$ and expensive physics predictions $\bm{c}_{2,s^{*}}$ are needed. We generate these predictions using the same Navier-Stokes-Boussinesq simulator but on coarser grids. More specifically, for each simulation, we run the fluid simulator from the same initialization as the ground truth but with a grid resolution of $32\times 32$. The buoyancy and vorticity fields at $s^{*}=10$ are the inexpensive physics prediction $\bm{c}_{1,s^{*}}$. Similarly, we run the fluid simulator with a grid resolution of $64\times 64$ and use the buoyancy as $\bm{c}_{2,s^{*}}$.

Four samples of $\bm{c}_{1,s^{*}}$ are plotted in the first column of Fig. 2. One can observe that the inexpensive simulation can capture the low-frequency patterns of the ground truth while details of fume swirls are blurred. This disparity demonstrates the inherent simulation bias.

In this study, we choose the conditional probability model $\log p(\bm{c}_{2,s^{*}}|\bm{x}_{0,s},\bm{c}_{1,s^{*}})$ as,

where $\text{AvgPool}^{2\times 2}$ is an average pooling operation (He et al., 2016) for the patch size of $2$ by $2$. More precisely, it divides the $64\times 64$ buoyancy field into $32\times 32$ patches of size $2\times 2$, then calculates the average buoyancy in each patch. Similarly, $\text{AvgPool}^{4\times 4}$ is a $4\times 4$ pooling operation. The model (5.1) encourages the low-frequency information of the simulated buoyancy and predicted buoyancy to be matched. $\gamma$ is a coefficient that measures the confidence of the result from the simulated buoyancy. We simply set it to be $0.01$ throughout our experiments. $C$ is a normalization constant that will become zero after differentiation. Additionally, we find removing the coefficient $1-\alpha_{t}$ in (Term 3) of (LABEL:eqn:expensivediffusionsample) and directly using $\bm{g}$ in the sampling update is more numerically stable. Therefore, we implement this modified version of the sampling update.

With access to physics simulations, we can implement MPDM to predict the ground-truth buoyancy. We use a U-Net (Ronneberger et al., 2015) implemented in (Wang, 2024) as the DeNN. The U-Net structure comprises eight 2D convolutional layers connected via skip connections. The diffusion step input $t$ first undergoes transformations through a sequence of sinusoidal functions with varying frequencies, followed by a nonlinear network to output temporal embeddings. These embeddings are then integrated into the intermediate layers of the U-Net through multiplication and addition operations. We use Adam to optimize the U-Net in Algorithm 2.

For comparisons, we also implement benchmark algorithms popular in the literature.

• •

  KOH (Kennedy & O’Hagan, 2000): The Bayesian calibration algorithm (KOH) uses a GP to model the residuals between the physics simulation output and the ground truth. We implement KOH to predict the difference between the ground truth buoyancy and the buoyancy from simulation using a batch-independent multi-output GP model provided by GPytorch (Gardner et al., 2018).

• •

  NN (Raissi et al., 2019): We directly train a deep neural network (a U-Net (Ronneberger et al., 2015)) to predict the future states of a system given historical information. The network is trained without inputs of $\bm{c}_{1,s^{*}}$ and $\bm{c}_{2,s^{*}}$.

• •

  Standard diffusion (S-DDIM): We implement the method described in Section 3.1 to sample target state vectors without the information from physics simulations.

The predictions of KOH, NN, the S-DDIM in Section 3.1, MPDM with inexpensive physics described in Section 3.2, and MPDM with expensive physics described in Section 3.3 are plotted in Fig. 2.

From Fig. 2, we can clearly see that S-DDIM predictions exhibit sharp and vivid details on the small-scale swirling structures of the fume. However, locations of large-scale swirl patterns are not accurate. This indicates that the DeNN learns the localized buoyancy patterns effectively but struggles to understand the long-range physical knowledge. The limitation is addressed when we use inexpensive physics prediction $\bm{c}_{1,s^{*}}$ as an additional input to the denoising network. Conditioned on the inexpensive physics prediction, samples from MPDM align more closely with the ground truth in terms of large-scale structures. The results imply that the physics knowledge is integrated in MPDM, which reaps benefits from the simulation while mitigating its bias. Furthermore, when $\bm{c}_{2,s^{*}}$ is available and used in choice 2 of Algorithm 3, the fidelity of the buoyancy prediction further improves, indicating that the more refined physics simulations can boost the quality of prediction.

On the contrary, the KOH model does not add meaningful information to the physics simulation, probably because the GP is known to have expressiveness limitations in the high-dimension regime (Wang et al., 2020). As our prediction is a $128\times 128=16382$ dimensional vector, predicting it using a Gaussian process is uneasy. The NN approach does not give accurate predictions either. This is understandable as NN completely relies on the training set to learn the physical evolution of the fume system. Such a purely statistical approach may not produce decent performance when the training set is not extremely large.

For numerical comparisons, we also evaluate the quality of the prediction on four standard evaluation metrics: mean squared error (MSE), peak signal-to-noise ratio (PSNR), structural similarity index (SSIM), and learned perceptual image patch similarity (LPIPS) (Zhang et al., 2018). In general, a lower MSE, a higher PSNR, a higher SSIM, and a lower LPIPS indicate a better sample quality. Amongst these metrics, MSE and PSNR measure the low-level pixel-wise difference between the sample and the ground truth, while SSIM is more aligned with human perception of the images. LPIPS leverages a deep neural network (VGG) to capture high-level visual similarities. The mean and standard deviation on the test set are reported in Table 1. In each column, the best result is highlighted in bold, while the second-best result is emphasized with an underline.

Results presented in Table 1 corroborate visual observations depicted in Fig. 2. S-DDIM is capable of generating samples that visually resemble the ground truth, as evidenced by its low LPIPS scores when compared with those from KOH and NN methods. However, the model’s high MSE and low PSNR indicate a poorer alignment with the ground truth at pixel levels.

With inexpensive physics predictions $\bm{c}_{1,s^{*}}$, MPDM enhances the sample quality as the LPIPS decreases. Furthermore, there are noticeable improvements in MSE, PSNR, and SSIM relative to the standard DDMs. Incorporation of the expensive physics simulation $\bm{c}_{2,s^{*}}$ into the reverse diffusion process (LABEL:eqn:expensivediffusionsample) leads to even greater improvements. The LPIPS scores decrease further, and the MSE, PSNR, and SSIM reach the highest values compared to all other evaluated algorithms. These improvements substantiate the benefits of utilizing multiple simulations in Algorithm 3.

In comparison, the benchmark KOH and NN incur high LPIPS, corroborating our observations that these samples do not show consistent visual patterns as the ground truth. The comparisons highlight the advantages of the proposed MPDM in producing high-quality results.

The melt pool often displays comet-like shapes, which is a result of heat dissipation and laser movement. In this section, we will construct a computationally amenable model to simulate the morphology and dynamics of the melt pool.

The interaction between high-energy laser beams and metal powders generates high-temperature particles that scatter from the metal surface. These particles can also be captured by the thermal camera. However, the 2D heat dissipation PDE can hardly describe the movement of particles. As a result, an alternative physics model is needed.

As the high-temperature particles are heated by the energy from the laser, it is natural to model them to be created at the center of the meltpool. After generation, these particles will move toward the edge of the receptive field of the thermal camera. The optical flow model (Horn & Schunck, 1981) provides a suitable characterization for such an emanating movement pattern.

More specifically, we define a velocity field $\bm{v}(x,y,s)=[v_{x},v_{y}]^{\top}\in\mathbb{R}^{2}$ denoting the velocity of the particle at position $(x,y)$ and time $s$. The velocity field $\bm{v}$ is the expensive physics simulation $\bm{c}_{2,s}$. We use $v_{x}$ to denote the velocity along the $x$-axis and $v_{y}$ to denote the velocity at the $y$-axis. For a small time interval $\Delta s$, the temperature of the scattering particles should remain approximately constant. As a result, $u_{p}$ should satisfy,

This equation is widely employed in the literature of optical flows (Horn & Schunck, 1981).

In practice, the optical flow may not be perfectly accurate because of the modeling, measurement, optimization, and statistical errors. We thus use the following probability model to characterize the inaccuracy,

where $\varepsilon$ are i.i.d. standard Gaussians noise. The parameter $\gamma$ represents our prior belief about the model accuracy: a larger $\gamma$ implies a higher confidence in the flow model. Therefore, given the optical flow $\bm{v}$, the conditional probability is,

Since we predict the temperature only on discrete pixels rather than the continuous 2D plane, the discrete counterpart of the equation above is,

where $\bm{v}$ plays the role of the expensive simulation $\bm{c}_{2,s^{*}}$.

In (6.2), $\mathcal{W}$ denotes the warping operator, which employs a semi-Lagrangian method to infer $\bm{x}_{0,s^{*}}$ from $\bm{x}_{0,s^{*}-1}$ and $\bm{v}$. A more detailed definition of $\mathcal{W}$ will be provided in the supplementary material. $\mathcal{P}_{\Omega_{\text{spatter}}}$ denotes the projection into the spatter region, $\Omega_{\text{spatter}}$, defined as follows,

Intuitively, $\mathcal{P}_{\Omega_{\text{spatter}}}$ applies a mask that selects the spatter temperature distribution from the entire temperature field. In our implementation, we define $\Omega_{\text{spatter}}$ as the complement of the melt pool region, which we estimate from the PDE solution (6.1.2).

It is important to highlight that both $\mathcal{W}$ and $\mathcal{P}_{\Omega_{\text{spatter}}}$ are differentiable. Hence, we can plug in (6.2) into (10) to calculate $\bm{g}$ using auto-differentiation packages. Subsequently, the gradient $\bm{g}$ is used in Algorithm 3 to guide the sampling process with the information from flow field $\bm{v}$.