Physics-Constrained Diffusion Model for Synthesis of 3D Turbulent Data

Abstract

Synthesizing fully developed three-dimensional turbulent velocity fields remains a long-standing problem in fluid mechanics and an open challenge for generative modeling. The difficulty arises from the coexistence of extreme dimensionality, multiscale rough fluctuations and strong intermittency, together with exact physical constraints such as incompressibility and zero-mean momentum. We propose a physics-constrained diffusion model (PCDM) in which these \emph{a priori} constraints are incorporated directly into the generative dynamics. Using rotating turbulence as a stringent benchmark, we show that the proposed framework enables stable and statistically faithful synthesis of inertial-range three-dimensional turbulent velocity fields at medium resolution, accurately reproducing anisotropic energy spectra, intermittency statistics, and physical constraints. By contrast, standard denoising diffusion probabilistic models without such constraints exhibit multiscale statistical deviations, violations of physical consistency, and substantially slower training convergence. These findings point to broader implications for generative modeling of high-dimensional complex systems under physical constraints.

Figures (9)

• Figure 1: Multiscale statistical structure of rotating turbulence. (a) Instantaneous three-dimensional rendering of the velocity magnitude $\|\boldsymbol{u}\|$, showing vertically aligned vortical structures induced by strong rotation. (b) Two-dimensional three-component (2D3C) component of the velocity field obtained from the decomposition of the full field in panel (a), dominated by coherent columnar structures. (c) Residual three-dimensional (3D) fluctuating field $\boldsymbol{u}_{3D}$ containing small-scale turbulent fluctuations; see Eq. (\ref{['eq:2d3c-3d']}) for details. (d) Total energy spectra for the full field (black solid line), the 2D3C component (red squares), and the 3D fluctuations (purple triangles). The shaded gray region indicates the forcing band around $k_f$, and the vertical dashed line marks the Kolmogorov wavenumber $k_\nu$. The inset shows the temporal evolution of the total kinetic energy normalized by its stationary mean; the shaded yellow interval denotes the portion of the simulation used to construct the training dataset (see Appendix \ref{['sec:rot_turb:dataset']} for details). (e) Compensated energy spectra $k^{5/3}E(k)$ highlighting the plateau associated with inertial-range scaling and the viscous cutoff (inset). (f) Probability density function (PDF) of the vertical vorticity component $\omega_z$, normalized by its standard deviation $\sigma(\omega_z)$, revealing pronounced non-Gaussian tails associated with intermittency. The inset shows a two-dimensional contour of $\omega_z$ in a plane perpendicular to the rotation axis, with red and blue indicating intense positive and negative vortical structures, respectively.
• Figure 2: Diffusion framework with physics-constrained backward dynamics (PCDM). (a) Forward and backward diffusion processes for a 3D turbulent velocity field. The forward process $q(\mathcal{V}_t|\mathcal{V}_{t-1})$ progressively injects Gaussian noise until the distribution approaches an isotropic Gaussian at $\mathcal{V}_T$, while the learned backward process $p_\theta(\mathcal{V}_{t-1}|\mathcal{V}_t)$ reverses this procedure to generate physically consistent velocity fields. (b) Physics-constrained backward transition in the PCDM. The neural network predicts the noise $\epsilon_\theta$, from which a clean estimate for the final generative step $t=0$, $\mathcal{V}_{0,\theta}$, is obtained. This estimated 3D field is projected onto the physically admissible subspace, enforcing incompressibility and zero-mean momentum, yielding $\tilde{\mathcal{V}}_{0,\theta}$. A sample $\mathcal{V}_{t-1}$ is then drawn from the Gaussian distribution defined by the reverse transition in Eq. (\ref{['eq:pcdm_reverse_process']}), and the procedure is repeated until the final sample $\mathcal{V}_0$ is obtained. The pseudocode below summarizes the training and sampling procedures of the PCDM. Removing the projection-based correction (highlighted in red in both the schematic and the algorithms) reduces the method to the standard DDPM formulation.
• Figure 3: Multiscale spectral diagnostics for DDPM-std, DDPM-prog, and the PCDM. Here DDPM-std and DDPM-prog denote the standard and progressively trained variants, respectively. (a) Energy spectra of the 2D3C component and the three-dimensional fluctuating field (3D FLUC), with the DNS reference (black solid line). (b) Ratio of generated to DNS spectra for the 2D3C component. (c) Same as (b), but for the 3D fluctuating component. Error bars are of the order of the symbol size and indicate variability estimated over three independent batches, each consisting of 200 generated fields.
• Figure 4: Training dynamics and convergence. (a) Pretraining loss of DDPM-prog, where the input domain is expanded from a thin slab toward progressively thicker domains. (b) Training loss of DDPM-std, DDPM-prog, and the PCDM on the full $64^3$ dataset. Three representative iterations, labeled A, B, and C, are selected from panel (b) to assess convergence. (c) Energy spectra of the 3D fluctuating component (3D FLUC) generated by DDPM-prog at iterations A, B, and C, compared with the DNS reference (black). (d) Same as (c), but for the PCDM. While both models exhibit progressive convergence during training, only the PCDM converges to the DNS spectrum across scales; DDPM-std exhibits persistent spectral fluctuations and is therefore not shown.
• Figure 5: Scale-dependent fourth-order flatness. Fourth-order flatness $K_\perp^{(4)}(r)$ of transverse velocity increments, shown as a function of the separation $r/l_\nu$, where $l_\nu$ denotes the Kolmogorov scale, measured in the plane normal to the rotation axis, for the reference DNS ($64^3$), DDPM-std, DDPM-prog, and the PCDM. Results are reported separately for (a) the slow 2D3C component and (b) the fast, fully three-dimensional fluctuating component (3D FLUC). The dashed horizontal line indicates the Gaussian reference value. The gray curve corresponds to the full-resolution DNS data ($256^3$; see Appendix \ref{['sec:rot_turb:dataset']}).