Paper deep dive
Neptuna: A Comprehensive Machine Learning Framework for Benchmarking Complex Multiphase Flows
Harish Ramachandran, Björn Kimpel, Thomas Paula, Josef Winter, Steffen Schmidt, Nikolaus Adams
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 89%
Last extracted: 8/2/2026, 1:21:23 PM
Summary
The paper introduces Neptuna, a comprehensive machine learning framework and benchmark for shock-driven compressible multiphase flows. It presents a 2.4 TB dataset of high-fidelity 2D and 3D simulations covering shock-induced bubble collapse and droplet breakup. The study evaluates diverse surrogate models (ConvNeXt, CNO, FFNO, ScOT, Poseidon, DPOT) using composite loss functions (MSE combined with H1-seminorm, SSIM, IRMSE) and adaptive weighting strategies (SoftAdapt, GradNorm). Results indicate that composite losses significantly improve interface preservation and spectral fidelity, with SoftAdapt providing the most consistent improvements.
Entities (15)
Relation Signals (11)
Neptuna → includes → 2.4 TB Dataset
confidence 95% · In this work, we introduce the first large-scale benchmark... comprising 2.4 TB of high-fidelity 2D and 3D datasets
SoftAdapt → improves → interface preservation
confidence 92% · SoftAdapt provides the most consistent improvements with almost no overhead compared to MSE-only training.
Neptuna → evaluates → ConvNeXt
confidence 90% · We evaluate diverse surrogate model families on our benchmarking framework: Neptuna... including ConvNeXt
Neptuna → evaluates → FFNO
confidence 90% · We evaluate diverse surrogate model families on our benchmarking framework: Neptuna... including... FFNO
Composite Losses → improves → spectral fidelity
confidence 90% · composite losses significantly improve interface preservation and spectral fidelity.
SABW → ispartof → 2.4 TB Dataset
confidence 90% · In this work, we introduce six datasets related to multiphase flows... Shock-induced Air Bubble Collapse in Water (SABW)
SDBA → ispartof → 2.4 TB Dataset
confidence 90% · Shock-induced Droplet Breakup in Air (SDBA)
SRBA → ispartof → 2.4 TB Dataset
confidence 90% · Shock-induced R22 Bubble Collapse in Air (SRBA)
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Compressible multiphase flows involving shocks and material interfaces arise in applications such as bubble collapse and droplet breakup, where strong nonlinear interactions produce complex interface deformation, mixing, and multiscale dynamics. Developing reliable machine learning surrogates for these flows remains challenging due to the simultaneous presence of compressibility, sharp discontinuities, and multiphase effects. In this work, we introduce the first large-scale benchmark specifically designed for shock-driven compressible multiphase flows, comprising 2.4 TB of high-fidelity 2D and 3D datasets featuring shock-induced bubble collapse and droplet breakup. We evaluate diverse surrogate model families on our benchmarking framework: Neptuna {this https URL}, including convolutional, spectral, transformer-based, and pre-trained PDE foundation models. Beyond standard MSE training, we investigate composite losses combining MSE with Sobolev, interface-aware, and structure-aware terms, together with adaptive loss balancing using SoftAdapt and GradNorm. Evaluation includes pointwise, spectral, feature-focused, structural, and physics-informed metrics. Results show that no single model performs best across all datasets and metrics, while composite losses significantly improve interface preservation and spectral fidelity. Among adaptive weighting strategies, SoftAdapt provides the most consistent improvements with almost no overhead compared to MSE-only training.
Tags
Links
- Source: https://arxiv.org/abs/2607.22280v2
- Canonical: https://arxiv.org/abs/2607.22280v2
Trouble viewing inline? Open PDF directly →
Full Text
177,493 characters extracted from source content.
Expand or collapse full text
Neptuna: A Comprehensive Machine Learning Framework for Benchmarking Complex Multiphase Flows Harish Ramachandran Technical University of Munich harish.ramachandran@tum.de Björn Kimpel Technical University of Munich bjoern.kimpel@tum.de Thomas Paula Technical University of Munich thomas.paula@tum.de Josef Winter Technical University of Munich josef.winter@tum.de Steffen Schmidt Technical University of Munich steffen.schmidt@tum.de Nikolaus Adams Technical University of Munich nikolaus.adams@tum.de Abstract Compressible multiphase flows involving shocks and material interfaces arise in applications such as bubble collapse and droplet breakup, where strong nonlinear interactions produce complex interface deformation, mixing, and multiscale dy- namics. Developing reliable machine learning surrogates for these flows remains challenging due to the simultaneous presence of compressibility, sharp disconti- nuities, and multiphase effects. In this work, we introduce the first large-scale benchmark specifically designed for shock-driven compressible multiphase flows, comprising 2.4 TB of high-fidelity 2D and 3D datasets 1 featuring shock-induced bubble collapse and droplet breakup. We evaluate diverse surrogate model families on our benchmarking framework: Neptuna 2 , including convolutional, spectral, transformer-based, and pre-trained PDE foundation models. Beyond standard MSE training, we investigate composite losses combining MSE with Sobolev, interface-aware, and structure-aware terms, together with adaptive loss balancing using SoftAdapt and GradNorm. Evaluation includes pointwise, spectral, feature- focused, structural, and physics-informed metrics. Results show that no single model performs best across all datasets and metrics, while composite losses sig- nificantly improve interface preservation and spectral fidelity. Among adaptive weighting strategies, SoftAdapt provides the most consistent improvements with almost no overhead compared to MSE-only training. 1 Introduction Compressible multiphase flows involving shock waves and material interfaces arise in a wide range of natural phenomena and engineering applications. Strong nonlinear interactions between shocks 1 Dataset repo: https://huggingface.co/FluidVerse. Dataset sample videos, metadata.json, inference rollout plots from autoregressive rollout of the trained baselines are provided in the supplementary_material.zip 2 Benchmarking repo: https://anonymous.4open.science/r/neptuna-A4E3 Preprint. arXiv:2607.22280v2 [physics.flu-dyn] 27 Jul 2026 and phase boundaries lead to complex dynamics such as interface deformation, vorticity generation, mixing, and topological changes. Shock-induced air bubble collapse in water represents an important regime of compressible multiphase flows. In gas–liquid configurations, such as air bubbles in mechanical equilibrium with water, the large acoustic impedance mismatch leads to strong compression and energy focusing during collapse, often producing high-pressure regions and jet formation. In contrast, gas–gas configurations, such as heavy gas (e.g. R22) bubbles in air, isolate compressibility-driven effects and provide a setup to analyze shock refraction, vorticity evolution, and mixing without the influence of surface tension or cavitation dynamics. Such phenomena have numerous real-world applications, including bubble curtains for protecting submerged structures (Würsig et al., 2000) and shockwave lithotripsy (Leighton, 2012). Another important class of problems is shock-induced droplet breakup in air. In this case, the interaction between the shock and a liquid droplet leads to a range of deformation and breakup mechanisms governed by the competition between inertial and surface tension forces. Depending on the Weber number, the droplet may undergo different breakup modes, such as Rayleigh Taylor Piercing (RTP) and Shear Induced Entrainment (SIE), each associated with a distinct instability mechanism (Theofanous and Li, 2008). These processes are highly relevant in applications such as fuel atomization, spray dynamics, and high-speed combustion (Lefebvre and McDonell, 2017). Despite significant progress in Computational Fluid Dynamics (CFD) solvers, simulating compress- ible multiphase flows with shocks remains computationally expensive due to the need for high spatial and temporal resolution to accurately capture the flow physics. These challenges have led to increas- ing interest in machine learning-based surrogate models for faster approximations. While several existing dataset repositories capture important aspects of compressible flow physics (Takamoto et al., 2022; Herde et al., 2024) or multiphase dynamics (Hassan et al., 2023, 2025), they do not jointly address compressibility, shock interactions, and multiphase effects in a comprehensive benchmark setting. To the best of our knowledge, this work provides the first benchmark datasets specifically designed for compressible multiphase flows with shocks, thereby extending the complexity of flow regimes in existing publically available datasets. We provide a 2.4 TB collection of high-fidelity 2D and 3D datasets and study reliable surrogate modeling for highly nonlinear, multiscale, and regime-dependent flows. We evaluate SOTA baselines across convolutional, spectral, and attention-based architectures, including ConvNeXt (CNeXt) (Liu et al., 2022), CNO (Raonic et al., 2023), FFNO (Tran et al., 2021), and ScOT (Herde et al., 2024), and fine-tune two pre-trained models: Poseidon (Herde et al., 2024) and DPOT (Hao et al., 2024). Beyond evaluating a broad spectrum of baselines, we train each model under two configurations: one using only Mean Squared Error (MSE), and another using a composite loss that augments MSE with additional terms designed to capture important flow characteristics. Specifically, we include gradient-based penalties via the H 1 -seminorm (Hartmann et al., 2019), interface-focused metrics such as the Interface-focused Root Mean Square Error (IRMSE) (Hassan et al., 2023), and structure-aware measures like the Structural Similarity Index (SSIM) (Nilsson and Akenine-Möller, 2020). These additional loss terms were selected empirically after extensive experimentation on the shock-induced air bubble collapse dataset, where their inclusion consistently improved either that said metric or a broad range of metrics during inference. To balance the contributions of these loss components, we employ two adaptive weighting strategies— SoftAdapt (SA) (Heydari et al., 2019) and Gradient Normalization (GN) (Chen et al., 2017)—which periodically sample each loss term and its gradient during training and adjust their respective weights accordingly. The framework also supports custom curriculum training blocks, enabling users to specify different training and evaluation strategies at different stages of training depending on the epoch. During inference, we report a comprehensive set of evaluation metrics beyond those used for training. These include additional pointwise measures such as Mean Absolute Error (MAE), Variance Scaled Root Mean Squared Error (VRMSE) (Ohana et al., 2024), and the H 2 -seminorm (Hartmann et al., 2019); spectral metrics like Multilevel Wavelet Loss (MLW) (Prantl et al., 2022) and wavelet frequency-binned RMSE (wfRMSE) (Takamoto et al., 2022); and feature-focused metrics that emphasize sharp gradients, such as shock-focused RMSE. We also include additional structure-aware measures like the Pearson Correlation Coefficient (PCC) (Hu et al., 2026), as well as physics-informed metrics that track integral quantities of interest, including mass, momentum, kinetic-energy for bubble datasets, center of mass displacement in case of droplets and enstrophy (Takamoto et al., 2022) over each rollout step. 2 A summarized contribution of this work is: •A 2.4 TB high-fidelity benchmark dataset for shock-driven compressible multiphase flows, covering bubble collapse and droplet breakup in 2D and 3D. •A systematic benchmark of diverse surrogate models, including convolutional, operator- based, transformer-based, and pre-trained PDE foundation models. •An investigation of physics- and structure-aware composite losses, combining MSE with derivative-, interface-, and structure-focused terms and adaptive loss balancing. • A reliability-focused evaluation framework extending beyond pointwise errors to spectral, structural, flow feature-focused, structure-aware and physics-informed metrics. 2 Related Work Developing reliable machine learning surrogates depends not only on the quality and diversity of the training data, but also on the model architecture, training strategy and loss functions used for optimization. Datasets capturing PDE dynamics continue to expand. Early publicly available collections, such as PDEBench (Takamoto et al., 2022) and PDEArena (Gupta and Brandstetter, 2022), focused on canonical PDEs that are relatively efficient to generate computationally. More recently, the field has shifted toward increasingly complex and sophisticated datasets that often require domain expertise for their curation, as showcased by the Well (Ohana et al., 2024), BlastNet (Chung et al., 2023), REALM (Mao et al., 2025), BubbleFormer (Hassan et al., 2025), and ClimSim (Yu et al., 2025). Additionally, RealPDEBench (Hu et al., 2026) presents experimental datasets aimed at fine-tuning pre-trained models, enabling them to better account for effects observed in real-world experiments. Many dataset papers now release their own benchmarking code alongside the data. Among recent efforts, APEBench (Koehler et al., 2024) provides more rigorous training strategies aimed at achieving stable long-horizon rollouts for neural operators. Related ideas also include perturbing inputs with Gaussian noise to improve robustness, as explored in Sanchez-Gonzalez et al. (2020). Several works have further introduced curriculum-style training strategies, such as the pushforward trick in Brandstetter et al. (2022), or training schedulers that transition efficiently between single-step and multistep autoregressive training over the course of epochs, as in Takamoto et al. (2023). Furthermore, there is a growing interest in training surrogate models with objectives that go beyond standard pointwise losses. Instead of relying only on MSE or MAE, several studies optimize weighted combinations of auxiliary, domain-specific metrics. For example, WaveLiT (Sankaran et al.) uses a combination of MSE and a wavelet-domain L 1 loss. PINO (Li et al., 2021) combines an L 2 loss term with a PDE-residual loss. Other works, such as Cho et al. (2024b), explore Sobolev losses during training, while CRONet (Olabiyi et al., 2025), from the structural topology optimization domain, uses SSIM loss to guide the optimization. 3 Datasets In this work, we introduce six datasets related to multiphase flows, generated by solving the compress- ible Euler equations. The generated datasets, as shown in Figure 1,span a wide variety of scenarios, including bubbles and droplets in both two and three dimensions, involving different materials and boundary conditions. To produce these datasets, we employ the high-fidelity finite volume solver ALPACA (Hoppe et al., 2022), incorporating the Robust Discrete Equations Method (RDEMIC) (Paula et al., 2023) for multiphase and interface treatment. Figure 2 shows the initial and boundary conditions used to generate the datasets with the dimensions mentioned in Table 1. The complete simulation domain is depicted by the black cuboid, while the orange cube denotes the cropped subdomain used for the dataset. We enrich the flow dynamics within the subdomain by considering multiple boundary-condition configurations. Symmetry boundaries behave as inviscid walls, reflecting the incident shock back into the subdomain. As an alternative, we obtain a subdomain free from boundary effects by running the simulation on a larger computational domain and subsequently cropping it to the desired region. We refer to these cropped, boundary- effect-free subdomains as having open boundaries. 3 Benchmarking Datasets AerobreakupShock Bubble Interaction Air Bubble in Water (SABW) Droplet in Air (SDBA) Varying Parameters Across Trajectories Droplet Count1 - 5 droplets Droplet Radius2m - 3.5m Droplet PositionRandom in domain Shock Mach Number1.2 - 3.2 Weber Number10 - 30 & 10k - 30k Varying Parameters Across Trajectories Bubble Count1 - 5 bubbles Bubble Radius2m - 3.5m Bubble PositionRandom in domain Shock Mach Number1.3 - 2.1 3D2D Boundary SSOOSS Boundary SSOO 128 3 1153 GB 256 2 52 GB 512 2 201 GB 144 train & 36 test trajectories each with 101 timesteps 240 train & 60 test trajectories each with 101 timesteps 144 train & 36 test trajectories each with 51 timesteps 2D Boundary O Boundary SSOOSS Boundary SSOO 3D 128 3 404 GB 256 2 40 GB 512 2 152 GB 256 2 39 GB 512 2 149 GB 240 train & 60 test trajectories each with 101 timesteps 240 train & 60 test trajectories each with 101 timesteps 240 train & 60 test trajectories each with 101 timesteps R22 Bubble in Air (SRBA) 2D Boundary O 256 2 43 GB 512 2 165 GB Figure 1: Dataset overview summarizing the provided datasets by dimensionality, boundary con- ditions, resolution, storage size, number of trajectories, and simulation-defined initial conditions. Boundary conditions are denoted using acronyms, whereSindicates symmetry andOindicates open boundaries. For example, in the 3D-SDBA dataset,SSOOSSrepresents symmetry boundaries on the North, South, Top, and Bottom faces, with open boundaries on the East and West. Y X Z Symmetry BC Zero- Gradient BC Symmetry BC Symmetry BC Symmetry BC Dirichlet BC Figure 2: Schematic of the dataset generation setup. The full simulation domain is outlined in black, while the orange cube highlights the extracted subdomain used as the provided dataset. The red plane indicates the planar shock wave which interacts with the downstream bubbles or droplets. 3.1 Shock-induced Air Bubble Collapse in Water (SABW) In the shock-induced air bubble collapse in water dataset, a moving planar shockwave impacts several bubbles resulting in complete breakdown of the bubbles into small fragments which get advected downstream and results in a series of interacting shockwaves emanating from the collapsed bubbles. A more detailed description of the physics of this dataset is provided in section A.3.1. The dataset is provided in two variants with different boundary conditions. In the first variant (O), all boundaries are set to open, enabling the study of shock–bubble interactions and inter-bubble dynamics in isolation. In the second variant (SSOO), symmetry boundary conditions are applied at the north and south boundaries, causing shockwave reflections, while the east and west boundaries remain open. This setup captures not only the direct interaction between the shock and the bubbles, as well 4 Table 1: Initial condition setup for different datasets with reference to Figure 2. Dataset L x L y L z S x d xw d xe d ys d yn d zt d zb [m][m][m][m][m][m][m][m][m][m] 2D-SABW O 76.896.0-25.828.210.228.828.8-- 2D-SABW SSOO 76.838.4-25.828.210.20.000.00-- 3D-SABW SSOOSS 76.838.438.425.828.210.20.000.000.000.00 2D-SRBA O 76.896.0-25.828.210.228.828.8-- 2D-SDBA SSOO 105.852.9-25.828.224.70.000.00-- 3D-SDBA SSOOSS 105.852.952.925.828.224.20.000.000.000.00 as inter-bubble effects, but also the influence of reflected shockwaves from the symmetric boundaries, acting as reflective walls. For the 3D dataset, only the symmetry-based variant (SSOOSS) is provided, as ensuring a sufficiently large simulation domain such that shockwaves from all boundaries do not reach the dataset region is computationally expensive. The temporal evolution of the density field for each variant is illustrated in Figures 3 and 4. Additional fields of interest, including pressure, velocity, volume fraction, and schlieren are provided in the Appendix (Figures 16 and 17). Density Figure 3: Uniformly spaced time snapshots of 2D-SABW at shock Mach number 1.70, with open boundary conditions on all sides (O), showing the evolution of density fromt 0 = 0.0 μsto t end = 15 μs. Density Figure 4: Uniformly spaced time snapshots of 2D-SABW at shock Mach number 1.70, with North–South symmetric and East–West open boundary conditions (SSOO), showing the evolu- tion of density from t 0 = 0.0 μs to t end = 15 μs. 3.2 Shock-induced Droplet Breakup in Air (SDBA) This dataset captures the dynamics of an external shock interacting with liquid droplets suspended in air. During the interaction, the shockwave is reflected, transmitted, and diffracted at the droplets, giving rise to surface instabilities. Depending on the Weber number—a dimensionless quantity representing the ratio of aerodynamic forces to surface tension, mainly determined by droplet diameter and surface tension coefficient—the interaction results in two distinct breakup regimes: Rayleigh–Taylor Piercing (RTP) (Figures 5 and 20), for Weber number below 100 and Shear-Induced Entrainment (SIE) (Figures 6 and 21) for Weber number above 1000. The detailed flow evolution and characterization of these regimes are provided in section A.3.2. Because droplet deformation requires long simulation times, implementing absorbing boundary conditions would require an excessively large computational domain to ensure that reflected shock waves do not re-enter the subdomain. Keeping computational constraints in mind, the dataset is generated using symmetry 5 boundary conditions, which are less computationally demanding, for both the two-dimensional SSOO and three-dimensional SSOOSS configurations. To further increase the complexity, both breakup modes are combined into a single dataset in equal proportion. This poses a significant challenge for the surrogate model, which must learn to infer the correct regime from the conditioning parameters during training and accordingly evolve its predictions at inference time. Density Figure 5: Uniformly spaced time snapshots of 2D-SDBA at shock Mach number 2.0, with all droplets in the RTP regime (Weber < 100) having a surface tension coefficientσ = 60.514, under North–South symmetric and East–West open boundary conditions (SSOO), showing the evolution of density from t 0 = 0.0 μs to t end = 215.5 μs. Density Figure 6: Uniformly spaced time snapshots of 2D-SDBA at shock Mach number 2.0, with all droplets in the SIE regime (Weber >1000) having a surface tension coefficientσ = 0.0719, under North–South symmetric and East–West open boundary conditions (SSOO), showing the evolution of density from t 0 = 0.0 μs to t end = 215.5 μs. 3.3 Shock-induced R22 Bubble Collaspe in Air (SRBA) In contrast to the SABW and SDBA datasets, the shock-interaction with R22 bubbles in air represents a compressible gas-gas configuration. This setup isolates density-driven interface dynamics and enables the study of shock refraction, baroclinic vorticity generation, and Richtmyer–Meshkov-type instabilities without the confounding effects of surface tension as observed in the droplet breakup datasets. Unlike the air bubble collapse, a mushroom shaped deformation is observed as seen in Figure 7 owing to the higher density of the R22 bubble. For this dataset, all boundaries are set to Open BC, a complete set of field variables is depicted in Figure 24 and the flow physics is described in section A.3.3. Density Figure 7: Uniformly spaced time snapshots of 2D-SRBA at shock Mach number 1.70, with open boundary conditions on all sides (O), showing the evolution of density fromt 0 = 0.0 μsto t end = 7.0 μs. 4 Benchmarking In this work, we provide a collection of trained two- and three-dimensional baseline models, including both randomly initialized networks and models fine-tuned from pre-trained weights. For the 2D setting, the randomly initialized baselines consist of 50M-parameter variants of ConvNeXt, CNO, 6 FFNO, and ScOT, while the pre-trained baselines include Poseidon-B (122M parameters) and DPOT- M (150M parameters). For the 3D setting, we benchmark 10M-parameter variants of FFNO and ConvNeXt. Both the 2D and 3D baselines are trained under a fixed 24-hour computational budget. Detailed model-specific hyperparameters are provided in Section B.1, while the overall training hyperparameters are summarized in Table 11. For models initialized from random weights, the learning problem is defined to extract four temporal snapshots from a dataset trajectory as input and predict the subsequent state. We further condition the model on the simulation parameters associated with the dataset, summarized in Figure 1, by using these parameters to modulate the normalization layers of the baseline architectures (Perez et al., 2017). The training and evaluation protocols are detailed in sections B.2 and B.3. For finetuning the pre-trained Poseidon-B model, we retain the original 1-to-1 input–output snapshot setting (Herde et al., 2024), while for finetuning DPOT-M, we follow the original 10-to-1 input–output configuration (Hao et al., 2024). 5 Results and Conclusion In this section, we present and compare benchmark results obtained using three training-loss strategies across the six datasets introduced in section 3. All baseline models described in section 4 are trained at a resolution of256 2 for the 2D datasets and128 3 for the 3D datasets. In our ablations, we consider three loss configurations: (i) training with mean squared error (MSE) alone; (i) training with a composite objective whose component weights are adapted using SA; and (i) training with the same composite objective, but with adaptive loss-weighting performed using GN. The composite objective combines MSE with additional terms designed to capture complementary aspects of the flow: the H 1 -seminorm (H1) penalizes errors in spatial gradients, the Structural Similarity Index (SSIM) encourages preservation of structure inside local windows and the interface-focused RMSE (IRMSE) emphasizes errors near bubble and droplet phase boundaries. Both SA and GN dynamically adjust the relative weights of these loss components during training. The additional computational overhead of these training strategies are analyzed from two perspectives: peak GPU memory requirement and time required for one epoch. The findings are reported in section C.1. We perform inference using a broader set of 12 evaluation metrics, as described in Table 13. Figures 8, 9, and 10 summarize representative results averaged across the test trajectories and the rollout steps for the three training-loss strategies. For each dataset, we report five of the twelve metrics, viz. VRMSE, MLW, IRMSE, SSIM and RMSE on the domain-integrated quantity of interest- enstrophy (iqRMSE-enstrophy). These metrics are selected to cover the five evaluation categories considered in this work: pointwise, spectral, feature-focused, structure-aware, and physics-informed metrics. To enable comparison across metrics with different scales, the raw metric values are first min-max normalized across all baseline models and training-loss strategies for a given dataset. The normalized values are then transformed by subtracting them from one, thereby reversing the scale so that lower original errors correspond to higher scores in the radar charts. Further details on the metrics are provided in section B.4, while the loss-weighting strategies are described in section B.5. The evolution of the errors as the baselines are autoregressively rolled out are presented in section C.2. The complete set of overall metrics aggregated across trajectories and rollouts, on all datasets and baselines are presented in section C.3. Selective baseline predictions during autoregressive rollout on the presented datasets are shown in section C.4. 5.1 Comparison of training loss strategies In this section we compare the results from the above-mentioned three different training loss strategies. Among all the 2D baselines, fine-tuned DPOT-M consistently outperforms the other baselines on the presented metrics. This could be attributed to either the rich pre-trained weights or the historic input sequence of 10 snapshots. We omit the fine-tuned Poseidon-B baseline from our discussion in this section because it performs substantially worse than the other baselines across nearly all metrics, despite having the largest parameter count. This could be attributed to the model not having a mechanism to incorporate conditioning, either through conditioning parameters or through inclusion of historic snapshots which help in distinguishing the different simulation trajectories. A further observation from Figures 8, 9, and 10 is that FFNO consistently outperforms CNeXt on the 2D bubble-collapse datasets. This suggests that spectral operator models may be better suited to 7 these cases, since the bubbles collapse within the first few timesteps and the dynamics are no longer dominated by a well-defined interface. As a result, the advantages of CNeXt as a convolutional edge detector become less pronounced. In contrast, for the shock-droplet breakup datasets, the droplets retain coherent interfaces as they deform and fragment, making the dynamics more localized and structure-preserving. This may explain why CNeXt performs slightly better than FFNO in these cases. Between the two composite-loss strategies, SA provides a more consistent improvement across the baselines and datasets compared to GN. VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SABW (O) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SABW (SSOO) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SRBA (O) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SDBA (SSOO) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 3D-SABW (SSOOSS) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 3D-SDBA (SSOOSS) MSE Training Loss ConvNeXtCNOFFNOScOTDPOT-MPoseidon-B Figure 8: Radar plot showcasing MSE-trained baselines across six datasets on five inference metrics- VRMSE, MLW, IRMSE, SSIM, and iqRMSE-enstrophy. VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SABW (O) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SABW (SSOO) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SRBA (O) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SDBA (SSOO) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 3D-SABW (SSOOSS) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 3D-SDBA (SSOOSS) Composite Training Loss with SoftAdapt Loss-Weighting Strategy ConvNeXtCNOFFNOScOTDPOT-MPoseidon-B Figure 9: Radar plot showcasing baselines trained with composite loss and SA loss-weighting across six datasets, evaluated on five inference metrics—VRMSE, MLW, IRMSE, SSIM, and iqRMSE- enstrophy VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SABW (O) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SABW (SSOO) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SRBA (O) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 2D-SDBA (SSOO) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 3D-SABW (SSOOSS) VRMSE MLW IRMSE SSIM iqRMSE-enstrophy 0.2 0.4 0.6 0.8 1.0 3D-SDBA (SSOOSS) Composite Training Loss with GradNorm Loss-Weighting Strategy ConvNeXtCNOFFNOScOTDPOT-MPoseidon-B Figure 10: Radar plot showcasing baselines trained with composite loss and GN loss-weighting across six datasets, evaluated on five inference metrics—VRMSE, MLW, IRMSE, SSIM, and iqRMSE- enstrophy. We now examine these trends in more detail by analyzing the representative metrics selected from each evaluation category for the 2D baselines. The models trained with MSE alone perform comparably to those trained with the two composite-loss strategies on the pointwise metric VRMSE (Equation 28). As shown in Tables 22, 27, 32, and 37, most baselines achieveVRMSEbelow 1. This indicates that the models are able to capture the mean background flow reasonably well, and that the composite objectives do not provide a substantial advantage for this pointwise metric. The effect of composite training is more pronounced for the interface-focused RMSE. For the bubble- collapse datasets, GN-based composite training improves IRMSE by approximately 5-7% on average across the baselines (Refer Tables 24, 29, 34). In contrast, for the shock-droplet interaction dataset, the improvement is substantially larger, reaching approximately 40-50% when using composite losses with SA and GN (Refer Table 39). This difference can be explained by the nature of the interface evolution in the two problem classes. After the bubble collapse, the interface, obtained by thresholding the density, becomes difficult to identify because the bubbles fragment into very small 8 structures. In the shock-droplet case, however, the droplet interface remains more coherent during deformation and breakup, making it easier to detect and therefore more responsive to interface-aware training. Significant improvements were observed for all baselines for the spectral metric-MLW across all baselines when trained with composite loss. There is also a small but consistent improvement in the high-frequency component of the wavelet frequency-binned RMSE when trained with composite loss as seen in Tables 23, 28, 33, 38. These improvements suggest that the metrics presented in our composite objective help reduce the low-frequency bias typically introduced by purely pointwise objectives, leading to better recovery of fine-scale structures. A consistent improvement is also observed for the structure-aware SSIM metric across all datasets and baselines, especially when using SA, followed by GN as shown in Tables 25, 30, 35, 40, 45 and 50. This trend may be partly attributed to the inclusion of SSIM as one of the terms in the composite training objective. However, we emphasize that incorporating a metric into the training objective does not necessarily guarantee improved performance on that metric at inference time; rather, the observed improvement in SSIM is an empirical outcome specific to the settings considered here. The improvement in SSIM metric during inference is also associated with better coherent structure prediction, as illustrated in Figure 36. Enstrophy of a flow field is computed as the integral summation of the absolute pointwise vorticities (Refer Table 19). For the domain integrated enstrophy metric (iqRMSE-enstrophy), we observe that there is no significant advantage of using composite losses while training. To a certain extent this holds true for other physics-based metrics presented in Tables 26, 31, 36, 41, 46 and 51. One exception is the ScOT baseline, which benefits substantially from composite training on selected datasets, particularly the 2D R22 bubble-collapse and 2D shock-droplet interaction cases. The 3D baselines generally perform worse than their 2D counterparts. As shown in the radar charts in Figures 8, 9, and 10, FFNO outperforms ConvNeXt across the reported metrics. Moreover, FFNO benefits consistently from composite-loss training, showing improvements across all selected metrics. However, these trends should be interpreted with caution. The 3D experiments were constrained by a 24-hour training limit, resulting in substantially fewer training epochs than in the 2D setting. More extensive training is required before drawing firm conclusions about the relative benefits of the composite losses in the 3D setting. 5.2 Conclusion In this work, we introduce six new datasets from the regime of compressible multiphase flows. We provide benchmarks on six baselines for the 2D datasets and two baselines for the 3D datasets. Each baseline is trained with three different loss strategies: MSE alone, and composite losses with adaptive loss-weighting using SoftAdapt and GradNorm. From the results, one can conclude that there is not a single baseline which performs best across all datasets and metrics and the choice of this surrogate would clearly depend on the metric of interest. However, the results do indicate that using composite losses with adaptive loss-weighting strategies can lead to substantial improvements across a range of metrics, especially those focused on interface quality and spectral content. Among the two adaptive loss-weighting strategies, SoftAdapt provides a more consistent improvement across the baselines and datasets compared to GradNorm, while also being more computationally efficient. In conclusion, our work highlights the importance of incorporating such complex flow datasets and diverse evaluation metrics that capture different aspects of the data into the SciML community to foster further advances in the development of robust data-driven surrogate models. Limitations and Future Work: There remains a significant gap in the ability of surrogate models to capture fine-scale details in complex flow problems. For example, none of the evaluated surrogate models were able to accurately recover the small droplets that detach from the surface of the primary droplet, although these secondary droplets are a defining characteristic of the SIE regime. Future work could therefore focus on improving the treatment of small-scale interface dynamics and fragmentation processes. Further improvements may also be achieved through more expressive conditioning strategies. In particular, incorporating information about boundary conditions during training could improve the generalizability of surrogate models across different simulation setups. The training metrics selected in this work for the composite losses were obtained after analyzing results from a broad list of experiments and selectively picking the metric to be included in the composite train 9 loss only when a positive correlation is obtained on the said metric or a diverse set of metrics during evaluation. In addition, a more detailed field-wise analysis of surrogate predictions would help identify specific failure modes and guide the development of improved architectures and training strategies. Finally, the benchmark could be extended to include more challenging datasets, such as reacting multiphase flows and multiphase flows with phase change, to further push the limits of surrogate modeling in this domain. 10 References John David Anderson. Modern compressible flow: with historical perspective. (No Title), 1990. Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton. Layer normalization. arXiv preprint arXiv:1607.06450, 2016. N Bempedelis and Y Ventikos. Energy focusing in shock-collapsed bubble arrays. Journal of Fluid Mechanics, 900:A44, 2020. Johannes Brandstetter, Daniel Worrall, and Max Welling. Message passing neural pde solvers. arXiv preprint arXiv:2202.03376, 2022. Zhao Chen, Vijay Badrinarayanan, Chen-Yu Lee, and Andrew Rabinovich. Gradnorm: Gradient normalization for adaptive loss balancing in deep multitask networks (2018). URL http://arxiv. org/abs/1711.02257. ArXiv, 1711, 2017. Namkyeong Cho, Junseung Ryu, and Hyung Ju Hwang. Sobolev training for operator learning, 2024a. URL https://arxiv.org/abs/2402.09084. Namkyeong Cho, Junseung Ryu, and Hyung Ju Hwang. Sobolev training for operator learning. arXiv preprint arXiv:2402.09084, 2024b. Wai Tong Chung, Bassem Akoush, Pushan Sharma, Alex Tamkin, Ki Sung Jung, Jacqueline Chen, Jack Guo, Davy Brouzet, Mohsen Talei, Bruno Savard, et al. Turbulence in focus: Benchmarking scaling behavior of 3d volumetric super-resolution with blastnet 2.0 data. Advances in Neural Information Processing Systems, 36: 77430–77484, 2023. Zachary Cooper-Baldock, Paulo E. Santos, Russell S. A. Brinkworth, and Karl Sammut. A generalised novel loss function for computational fluid dynamics, 2024. URL https://arxiv.org/abs/2411.17059. Nico Fleischmann, Stefan Adami, and Nikolaus A Adams. A shock-stable modification of the hllc riemann solver with reduced numerical dissipation. Journal of computational physics, 423:109762, 2020. Sigal Gottlieb and Chi-Wang Shu. Total variation diminishing runge-kutta schemes. Mathematics of computation, 67(221):73–85, 1998. Jayesh K Gupta and Johannes Brandstetter. Towards multi-spatiotemporal-scale generalized pde modeling. arXiv preprint arXiv:2209.15616, 2022. Zhongkai Hao, Chang Su, Songming Liu, Julius Berner, Chengyang Ying, Hang Su, Anima Anandkumar, Jian Song, and Jun Zhu. Dpot: Auto-regressive denoising operator transformer for large-scale pde pre-training. arXiv preprint arXiv:2403.03542, 2024. Thomas Hartmann, Matthias Kissel, and Klaus Diepold. Sobolev training with higher order derivatives. Technical report, Chair of Data Processing, Technical University of Munich, 2019. URLhttps://collab.dvb.bayern/ download/attachments/75112352/fp_hartmann.pdf. Accessed: 2026-03-05. Sheikh Md Shakeel Hassan, Arthur Feeney, Akash Dhruv, Jihoon Kim, Youngjoon Suh, Jaiyoung Ryu, Yoonjin Won, and Aparna Chandramowlishwaran. Bubbleml: A multi-physics dataset and benchmarks for machine learning. arXiv preprint arXiv:2307.14623, 2023. Sheikh Md Shakeel Hassan, Xianwei Zou, Akash Dhruv, and Aparna Chandramowlishwaran. Bubbleformer: Forecasting boiling with transformers. In The Thirty-ninth Annual Conference on Neural Information Processing Systems Datasets and Benchmarks Track, 2025. Maximilian Herde, Bogdan Raonic, Tobias Rohner, Roger Käppeli, Roberto Molinaro, Emmanuel de Bézenac, and Siddhartha Mishra. Poseidon: Efficient foundation models for pdes. Advances in Neural Information Processing Systems, 37:72525–72624, 2024. A Ali Heydari, Craig A Thompson, and Asif Mehmood. Softadapt: Techniques for adaptive loss weighting of neural networks with multi-part loss functions. arXiv preprint arXiv:1912.12355, 2019. Nils Hoppe, Josef M Winter, Stefan Adami, and Nikolaus A Adams. Alpaca-a level-set based sharp-interface multiresolution solver for conservation laws. Computer Physics Communications, 272:108246, 2022. Peiyan Hu, Haodong Feng, Hongyuan Liu, Tongtong Yan, Wenhao Deng, Tianrun Gao, Rong Zheng, Haoren Zheng, Chenglei Yu, Chuanrui Wang, et al. Realpdebench: A benchmark for complex physical systems with real-world data. arXiv preprint arXiv:2601.01829, 2026. 11 Guang-Shan Jiang and Chi-Wang Shu. Efficient implementation of weighted eno schemes. Journal of computa- tional physics, 126(1):202–228, 1996. Felix Koehler, Simon Niedermayr, Rüdiger Westermann, and Nils Thuerey. Apebench: A benchmark for autoregressive neural emulators of pdes. Advances in Neural Information Processing Systems, 37:120252– 120310, 2024. Olivier Le Métayer and Richard Saurel. The noble-abel stiffened-gas equation of state. Physics of Fluids, 28(4), 2016. Arthur H Lefebvre and Vincent G McDonell. Atomization and sprays. CRC press, 2017. Timothy Leighton. The acoustic bubble. Academic press, 2012. Z Li, H Zheng, N Kovachki, D Jin, H Chen, B Liu, K Azizzadenesheli, and A Anandkumar. Physics-informed neural operator for learning partial differential equations. arxiv. arXiv preprint arXiv:2111.03794, 2021. Ze Liu, Han Hu, Yutong Lin, Zhuliang Yao, Zhenda Xie, Yixuan Wei, Jia Ning, Yue Cao, Zheng Zhang, Li Dong, et al. Swin transformer v2: Scaling up capacity and resolution. 2022 ieee. In CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 11999–12009, 2021. Zhuang Liu, Hanzi Mao, Chao-Yuan Wu, Christoph Feichtenhofer, Trevor Darrell, and Saining Xie. A convnet for the 2020s. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 11976–11986, 2022. Runze Mao, Rui Zhang, Xuan Bai, Tianhao Wu, Teng Zhang, Zhenyi Chen, Minqi Lin, Bocheng Zeng, Yangchen Xu, Yingxuan Xiang, et al. Benchmarking neural surrogates on realistic spatiotemporal multiphysics flows. arXiv preprint arXiv:2512.18595, 2025. Jomela C Meng and Tim Colonius. Numerical simulation of the aerobreakup of a water droplet. Journal of Fluid Mechanics, 835:1108–1135, 2018. Jim Nilsson and Tomas Akenine-Möller. Understanding ssim. arXiv preprint arXiv:2006.13846, 2020. Robert R Nourgaliev, Truc-Nam Dinh, and Theo G Theofanous. Adaptive characteristics-based matching for compressible multifluid dynamics. Journal of Computational Physics, 213(2):500–529, 2006. Ruben Ohana, Michael McCabe, Lucas Meyer, Rudy Morel, Fruzsina J Agocs, Miguel Beneitez, Marsha Berger, Blakesley Burkhart, Stuart B Dalziel, Drummond B Fielding, et al. The well: a large-scale collection of diverse physics simulations for machine learning. Advances in Neural Information Processing Systems, 37: 44989–45037, 2024. Ridwan Olabiyi, Hui Yang, and Ashif Iquebal. Cronet: A convolutional recurrent operator approximator network to accelerate topology optimization. Manufacturing Letters, 44:1052–1063, 2025. Thomas Paula, Stefan Adami, and Nikolaus A Adams. A robust high-resolution discrete-equations method for compressible multi-phase flow with accurate interface capturing. Journal of Computational Physics, 491: 112371, 2023. E Perez, F Strub, H De Vries, and V Dumoulin. Visual reasoning with a general conditioning layer, courville. In In Proceedings of the AAAI Conference on Artificial Intelligence, 2017. Lukas Prantl, Jan Bender, Tassilo Kugelstadt, and Nils Thuerey. Wavelet-based loss for high-frequency interface dynamics, 2022. URL https://arxiv.org/abs/2209.02316. Bogdan Raonic, Roberto Molinaro, Tobias Rohner, Siddhartha Mishra, and Emmanuel de Bezenac. Convolutional neural operators. In ICLR 2023 workshop on physics for machine learning, 2023. Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International conference on machine learning, pages 8459–8468. PMLR, 2020. Shyam Sankaran, Hanwen Wang, and Paris Perdikaris. Wavelit: A parameter-efficient architecture for neural pde solvers. Makoto Takamoto, Timothy Praditia, Raphael Leiteritz, Daniel MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert. Pdebench: An extensive benchmark for scientific machine learning. Advances in neural information processing systems, 35:1596–1611, 2022. 12 Makoto Takamoto, Francesco Alesiani, and Mathias Niepert. Learning neural pde solvers with parameter-guided channel attention. In International Conference on Machine Learning, pages 33448–33467. PMLR, 2023. TG Theofanous and GJ Li. On the physics of aerobreakup. Physics of fluids, 20(5), 2008. Eleuterio F Toro. Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science & Business Media, 2013. Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized fourier neural operators. arXiv preprint arXiv:2111.13802, 2021. Abhinau K Venkataramanan, Chengyang Wu, Alan C Bovik, Ioannis Katsavounidis, and Zafar Shahid. A hitchhiker’s guide to structural similarity. IEEE Access, 9:28872–28896, 2021. Josef Winter, Jakob Kaiser, Stefan Adami, and Nikolaus Adams. Numerical investigation of 3d drop-breakup mechanisms using a sharp interface level-set method. In 11th International Symposium on Turbulence and Shear Flow Phenomena, TSFP 2019, 2019. B Würsig, CR Greene Jr, and TA Jefferson. Development of an air bubble curtain to reduce underwater noise of percussive piling. Marine environmental research, 49(1):79–93, 2000. Sungduk Yu, Zeyuan Hu, Akshay Subramaniam, Walter Hannah, Liran Peng, Jerry Lin, Mohamed Aziz Bhouri, Ritwik Gupta, Björn Lütjens, Justus C Will, et al. Climsim-online: A large multi-scale dataset and framework for hybrid physics-ml climate emulation. Journal of Machine Learning Research, 26(142):1–85, 2025. 13 Appendix Table of Contents A Dataset Details15 A.1 Dataset Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .15 A.2 Dataset Specifications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .15 A.3 Base Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .16 A.4 Downsampling of Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . .31 B Benchmarking Details36 B.1 Baseline Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .36 B.2 Training Protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .38 B.3 Evaluation Protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .39 B.4 Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .39 B.5 Adaptive loss balancing framework . . . . . . . . . . . . . . . . . . . . . . . .47 B.6 Computational Resources . . . . . . . . . . . . . . . . . . . . . . . . . . . . .49 C Detailed results50 C.1 Metric performance comparison . . . . . . . . . . . . . . . . . . . . . . . . . .50 C.2 Inference metric evolution over rollout steps . . . . . . . . . . . . . . . . . . .50 C.3 Overall inference metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . .76 C.4 Rollout visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .89 14 A Dataset Details The datasets are made publicly available on HuggingFace athttps://huggingface.co/FluidVerse. A.1 Dataset Organization Each dataset is divided into training and testing sets and stored in HDF5 format using a hierarchical organization. At the top level, groups are defined, with their names encoding the parameters that specify the initial conditions of each trajectory. Within each group, the time-series field data is stored in the shapeT × C × X res × Y res × Z res , whereT, C, X res , Y res , Z res denote the number of timesteps, channels, and resolution in the X, Y and Z direction respectively. An illustration of this hierarchical structure is provided below: Group: 005_Mas1.30_sb1_Ax0.0078_Ay0.0179_Az0.0145_Ar0.0022 Field: density Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: density_1 Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: density_2 Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: diffuse_volume _fraction_1 Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: pressure Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: velocityX Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: velocityY Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: velocityZ Shape: (51, 1, 128, 128, 128) Dtype: float64 . . . Group: 009_Mas1.30_sb1_Ax0.0130_Ay0.0309_Az0.0097_Ar0.0025 Field: density Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: density_1 Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: density_2 Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: diffuse_volume _fraction_1 Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: pressure Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: velocityX Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: velocityY Shape: (51, 1, 128, 128, 128) Dtype: float64 Field: velocityZ Shape: (51, 1, 128, 128, 128) Dtype: float64 During dataloading, the simulation parameters from the group name are extracted, scaled using min- max normalization and passed to the normalization layer of the baseline architecture as conditioning parameters as described in section B.2. A.2 Dataset Specifications As mentioned in the main text, the datasets are generated using the high-fidelity finite volume solver ALPACA. Within this framework, the Robust Discrete Equations Method for Interface Cap- turing (RDEMIC) is employed to produce ground-truth trajectories by solving the two- and three- dimensional compressible Euler equations. The governing equations in vector notation are given as follows: ∂ t U l +∇· F l = B l ·∇α l (1) where subscriptldenotes the index of the phase,U l is the vector of conserved quantities,F l is the flux tensor and B l is the interaction tensor. 15 U l = α l α l ρ l α l ρ l u l α l E l ,F l = 0 α l ρ l u T l α l ρ l u l ⊗ u l + α l p l I α l (E l + p l )u T l ,B l = −u T int 0 p int,l I p int,l u T int , The quantitiesα l ,ρ l ,u l ,p l , andE l represent the volume fraction, density, velocity vector, pressure, and total energy of phasel, respectively. The interface velocity and pressure are denoted byu int and p int,l . In the absence of surface tension, the interface pressurep int,l is identical across all phases.I represents the identity tensor. To close the above governing equations, an additional Equation of State (EOS) is required which relates pressure, density, and internal energy. In this work, the stiffened-gas EOS is employed for the generation of all datasets, and is given by p(ρ,e) = (γ− 1)ρe− γp stiff (2) wherepdenotes the fluid pressure,ρthe mass density,ethe specific internal energy, andγa model constant. The parameter p stiff represents the contribution from fluid pre-compression. The simulations use the HLLC Riemann solver (Toro, 2013). In case of air-bubble collapse in water, the HLLC-LM solver is used, which remedies shock instabilities associated with the low- Mach-number flow in transverse direction to the shock wave (Fleischmann et al., 2020). A WENO5 spatial reconstruction (Jiang and Shu, 1996) is employed to achieve high accuracy and robust shock capturing. Time integration is carried out using a third-order Runge-Kutta scheme (Gottlieb and Shu, 1998), ensuring stable and accurate evolution of the flow physics. Table 2 presents the summarized metadata for the generated datasets. The resolution listed here refers to the original output generated directly by the solver; the data is subsequently cropped and downsampled to the target resolution specified in Figure 1. Table 2: Metadata table for the provided datasets. End-time refers to the total simulation time, CFL refers to Courant-Friedrichs-Lewy criterion,∆t solver is the average solver timestep across trajectories and∆xrefers to the grid spacing. The reported wall-clock time is obtained using the compute resources described in section B.6. DatasetResolution End-time [μs] CFL ∆t solver [μs] ∆x [m] # snap- shots Wall-clock time [hours] 2D-SABW O [2048× 2560]150.51.6×10 −3 0.037510146 2D-SABW SSOO [2048× 1024]150.51.3×10 −3 0.037510127 3D-SABW SSOOSS [512× 256× 256]150.54.35×10 −3 0.1500511400 2D-SRBA O [2048× 2560]700.513.0×10 −3 0.037510121 2D-SDBA SSOO [1024× 512]91 - 12400.413.0×10 −3 0.103010146 3D-SDBA SSOOSS [256× 128× 128]91 - 12400.436.8×10 −3 0.413351620 A.3 Base Setup The datasets capture the time-evolving behavior of bubbles and droplets subjected to an external shock wave. The computational domain is partitioned into pre-shock and post-shock regions. Once the shock wave interacts with the bubble or droplet, the resulting post-shock flow largely governs their deformation and subsequent dynamics. The shock Mach number (M s ) serves as an influential non-dimensional parameter characterizing this post-shock regime. For each initial condition, the 16 post-shock density, velocity and pressure are set as Dirichlet boundary conditions in the left boundary of the domain. A.3.1 Shock-induced Air Bubble Collapse in Water (SABW) In a multiphase flow, one of the key important fields is the evolution of the density. Here, we examine its evolution when an external shock wave interacts with a cylindrical air bubble in water, as shown in Figure 11. Att = 1.5 μs, the shock wave reaches the bubble, causing the air inside to accelerate due to the transmitted shock. Subsequently, a water jet begins to penetrate the bubble from the upstream side, as observed att = 3.9 μs. Byt = 4.8 μs, this jet impacts the downstream interface of the bubble, generating a localized peak pressure in the domain, commonly referred to as the water hammer effect. This phenomenon produces a radially propagating shock wave and leads to the breakup of the bubble into two fragments in the two-dimensional cylindrical case, or into a toroidal structure in the three-dimensional spherical case. The resulting fragments undergo further breakup under the influence of the radial shock, producing additional pressure waves, as seen att = 8.4 μs. Eventually, the fragmented bubble is advected downstream with the flow. Figure 11: Time evolution of density field for the 2D-SBAW-O dataset. The post-shock properties for the SABW configuration are determined using the normal shock relations (Le Métayer and Saurel, 2016), as given in Equation 3. As discussed in the base setup, the domain is divided into pre- and post-shock regions; in the case of SABW, the fluid in the domain is water, and thus the properties of water are used to compute the pre-shock values. In this formulation, M s denotes the shock Mach number, whilec pre represents the speed of sound in the pre-shock region, calculated using the stiffened gas equation. The pre-shock density (ρ pre ) and pressure (p pre ) are specified as 1000 kg/m 3 and 1.0 × 10 5 Pa, respectively. Additionally, for water, the background pressure (p ∞ ) and the ratio of specific heats (γ) are taken as 60 × 10 5 Pa and 4.4. The density of the bubble is initialized to 1.0 kg/m 3 and the pressure inside is initialized to 1.0 × 10 5 Pa. 17 c pre = s γ p pre + p ∞ ρ pre u s = M s ∗ c pre p post = (p pre + p ∞ ) 1 + 2γ (γ + 1) (M 2 s − 1) − p ∞ ρ post = ρ pre (γ + 1)M 2 s 2 + (γ− 1)M 2 s u post = u s 1− ρ pre ρ post (3) The solver was validated on two setups, a single bubble setup and a multi bubble setup. For the single bubble setup, we validate the evolution of the rate of change of volume of air over a given time interval and compare it against the results shown in Nourgaliev et al. (2006). The setup and the corresponding geometric parameters are shown in Figure 12 and Table 3 respectively. The shock Mach number is set to 1.72, the CFL number is 0.4 and the bubble and the surrounding are initialized at atmospheric pressure of 1 bar. The simulation is run for a total duration of 5μs. We perform the simulation at two different resolutions of512× 512and256× 256such that there are 100 cells and 50 cells, respectively, in the inital diamter of the bubble. Figure 13 shows the rate of change of volume over time for both resolutions and the results are in good agreement with the reference plot sourced from Nourgaliev et al. (2006). Shock wave Air Post- Shock Region (Water) Pre-Shock Region (Water) Symmetry BC Symmetry BC Zero-Gradient BC Zero-Gradient BC Figure 12: Single air bubble in water setup. The parameter details are presented in Table 3. The shock Mach number is 1.72 and the bubble and the surrounding are initialized at atmospheric pressure of 1 bar. Table 3: Geometric parameters for the single air bubble in water setup. ParameterValue [m] L x 0.03072 L y 0.03072 S x 0.01040 r c 0.00300 (x c ,y c )(0.01548, 0.001536) The next validation case evaluates the solver by comparing the peak pressure in a multi-bubble configuration with a well-established reference (Bempedelis and Ventikos, 2020). As illustrated in 18 012345 Time [ μs ] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 dV dt [ m 3 s ] Reference Data 50 Cells/Dia Uniform Grid 100 Cells/Dia Uniform Grid Figure 13: Comparing evolution of rate of change of volume with time for the single air bubble in water with the reference plot shown in Nourgaliev et al. (2006). Figure 14, the setup consists of three staggered air bubbles immersed in water, with the corresponding geometric parameters listed in Table 4. The shock Mach number is set to 1.42, the CFL number to 0.5, and both the bubbles and the surrounding medium are initialized at an atmospheric pressure of 1 bar. The simulation is carried out over a total duration of 1.5μs. Validation is achieved by comparing the temporal evolution of the peak pressure in the domain with the results reported in Bempedelis and Ventikos (2020). The comparison shows good agreement with the reference, with the solver accurately capturing the key events over time. The slight discrepancy observed in the maximum peak pressure throughout the simulation is attributed to differences in interface treatment: the reference employs a level-set method, whereas the present solver uses a diffuse interface approach, which leads to some smearing of peak values. Periodic BC Periodic BC Zero-Gradient BC Zero-Gradient BC Post- Shock Region (Water) Pre-Shock Region (Water) Shock wave Air Air Air Figure 14: Three staggered air bubble in water setup. The parameter details are presented in Table 4. The shock Mach number is set to 1.42 and the bubbles and the surrounding are initialized at atmospheric pressure of 1 bar. 19 Table 4: Geometric parameters for the staggered three air bubbles in water setup. ParameterValue [m] L x 0.00350 L y 0.00200 S x 0.00015 R b 0.00050 r b 0.00020 (d x ,d y )(0.00050,0.00150) 0.00.20.40.60.81.01.21.41.6 Time [ μs ] 0.0 0.5 1.0 1.5 2.0 2.5 Pressure [ Pa ] ×10 10 Reference Simulation Figure 15: Peak pressure variation with time for the staggered three air bubble in water. The resolution for the simulation is 896× 512. 20 Density Pressure Velocity X Velocity Y Volume Fraction Schlieren Figure 16: Uniformly spaced time snapshots of Shock-induced Air Bubble collapse in Water (SABW) at shock Mach number 1.70, with open boundary conditions on all sides (O), showing the evolution of pressure, density, x-velocity, y-velocity, volume fraction, and schlieren fromt 0 = 0.0 μs to t end = 15 μs. 21 Density Pressure Velocity X Velocity Y Volume Fraction Schlieren Figure 17: Uniformly spaced time snapshots of Shock-induced Air Bubble collapse in Water (SABW) at shock Mach number 1.70, with North–South symmetric and East–West open boundary conditions (SSOO), showing the evolution of pressure, density, x-velocity, y-velocity, volume fraction, and schlieren from t 0 = 0.0 μs to t end = 15 μs. 22 |velocity| vorticity pressure schlieren Figure 18: Uniformly spaced time snapshots of a 3D shock-induced air bubble collapse in water (SABW) at a shock Mach number of 1.70. The configuration employs symmetry boundary conditions in the North–South and Top–Bottom directions and open boundaries in the East–West direction (SSOOSS). The figure illustrates the temporal evolution from t 0 = 0.0 μs to t end = 15 μs, showing pressure–vorticity fields on the horizontal plane, schlieren–velocity magnitude on the vertical plane, and a yellow density iso-contour of the bubbles. 23 A.3.2 Shock-induced Droplet Breakup in Air (SDBA) The shock induced droplet breakup datasets constitute two distinct breakup modes, namely Rayleigh Taylor Piercing (RTP) and Shear Induced Entrainmet (SIE). The transition between these modes is primarily determined by the Weber number (We), which represents the ratio of aerodynamic forces to surface tension. We = ρ post u 2 post D σ (4) In the above equation, D is the droplet diameter and σ is the surface tension coefficient. When a shock wave impacts a droplet, the initial response—largely independent of the Weber number—is a deformation phase in which the droplet flattens. This deformation arises from a non- uniform pressure distribution around its surface. The subsequent stage is the actual breakup, which is strongly governed by the Weber number. The RTP breakup regime typically occurs for Weber numbers between roughly 10 and 100, whereas the SIE regime dominates at much higher values (We > 1000), with intermediate Weber numbers representing a transition between the two modes (Theofanous and Li, 2008). In the SIE regime, breakup is driven mainly by strong shear forces acting along the droplet surface. After the droplet has flattened, shear-induced disturbances emerge near the equator; these instabilities originate near the droplet equator after the droplet has flattened out in the first phase and are advected along the droplet surface. As the relative velocity between the droplet and the surrounding gas increases, these disturbances grow due to Kelvin–Helmholtz instability, eventually stripping liquid from the droplet and producing fine droplets downstream. In contrast, in the RTP regime, relatively stronger surface tension suppresses the growth of such shear instabilities, maintaining a smoother interface. As deformation progresses, the upstream side of the droplet becomes concave as the surrounding gas penetrates and pierces the liquid. Unlike RTP, the SIE regime is characterized by a continuous and gradual loss of mass, often resulting in a mist of droplets downstream (Winter et al., 2019). For determining the post-shock properties for the SDBA setup, we use the shock relations as shown in Equation 3 to compute the flow properties behind the shock (Anderson, 1990). u s = M s · c pre u pre,rel =−u s u pre = u pre,rel + u s T post = T pre 1 + 2γ M 2 s − 1 γ + 1 ! 2 + (γ− 1)M 2 s (γ + 1)M 2 s c post = p γ· R· T post M f,post,rel = s 1 + γ−1 2 M 2 s γM 2 s − γ−1 2 u post,rel = M f,post,rel · c post u post = u s − u post,rel ρ post = ρ pre · (γ + 1)M 2 s 2 + (γ− 1)M 2 s p post = p pre 1 + 2γ M 2 s − 1 γ + 1 ! (5) In the above equation set,M f is the post-shock flow Mach number,Tis the temperature,cis the speed of sound,γ = c p c v is the ratio of specific heat capacities set to 1.4, andRis the specific gas 24 constant, set to 287 J/kg K. The pre-shock density,ρ pre , = 1.2 kg/m 3 , the pre-shock pressure,p pre , = 101325 Pa and the pre-shock temperature,T pre = 300 K. The density of the droplet is 1000 kg/m 3 and pressure initialized inside the droplet is 101325 Pa. The solver was validated for the droplet scenario by comparing the center of mass displacement with the expected trend as shown in Meng and Colonius (2018) when a 2D-axissymmetric spherical droplet whose center is placed on the x-axis interacts with a shockwave of Mach number 1.20. 0.000.250.500.751.001.251.501.752.00 Time ×10 3 [ s] 0.040 0.045 0.050 0.055 0.060 0.065 0.070 0.075 x-COM displacement [m] Droplet x-COM displacement vs Time Weber number - 38939, Droplet diamter- 4m, Mach number - 1.20 Simulation data Expected trend Figure 19: Center-of-mass displacement in x-axis when a 4-m droplet interacts with a shockwave of Mach 1.2 A.3.3 Shock-induced R22 Bubble Breakup in Air (SRBA) The dynamics for the SRBA case is notably different from the SABW scenario. Since the R22 bubble is denser than the surrounding air, the transmitted shock within the bubble propagates more slowly than the incident shock, as illustrated in Figure 23(a). As the transmitted shock converges toward the downstream side of the bubble, a small cusp forms, as seen in Figure 23(b), which eventually evolves into a thin R22 jet, shown in Figure 23(c). During the same time interval, small vortices develop along the interface and are subsequently amplified as they move downstream of the bubble. For determining the post-shock properties for SRBA setup, we again use the normal shock relations given in Equation 3. In this case, the medium in the domain is air, and thus the properties of air are used to compute the pre-shock values. The pre-shock density (ρ pre ) and pressure (p pre ) are specified as 1.2041 kg/m 3 and 1.01×10 5 Pa, respectively. Furthermore, for air, the background pressure (p ∞ ) and the ratio of specific heats (γ) are taken as 0.00 Pa and 1.249. The R22 bubble is initialized with a density of 3.7975 kg/m 3 . 25 Density Pressure Velocity X Velocity Y Volume Fraction Schlieren Figure 20: Uniformly spaced time snapshots of Shock-induced Droplet Breakup in Air (SDBA) at shock Mach number 2.0, all droplets in the RTP regime (Weber < 100) having a surface tension coefficientσ = 60.514, under North–South symmetric and East–West open boundary conditions (SSOO), showing the evolution of pressure, density, x-velocity, y-velocity, volume fraction, and schlieren from t 0 = 0.0 μs to t end = 215.5 μs. 26 Density Pressure Velocity X Velocity Y Volume Fraction Schlieren Figure 21: Uniformly spaced time snapshots of Shock-induced Droplet Breakup in Air (SDBA) at shock Mach number 2.0, with all droplets in the SIE regime (Weber > 1000) having a surface tension coefficientσ = 0.0719, under North–South symmetric and East–West open boundary conditions (SSOO), showing the evolution of pressure, density, x-velocity, y-velocity, volume fraction, and schlieren from t 0 = 0.0 μs to t end = 215.5 μs. 27 |velocity| vorticity pressure schlieren Figure 22: Uniformly spaced time snapshots of a 3D shock-induced droplet breakup in air (SDBA) at a shock Mach number of 1.20, with the breakup mode in the SIE regime. The configuration employs symmetry boundary conditions in the North–South and Top–Bottom directions and open boundaries in the East–West direction (SSOOSS). The figure illustrates the temporal evolution fromt 0 = 0.0μs tot end = 1240 μs, showing pressure–vorticity fields on the horizontal plane, schlieren–velocity magnitude on the vertical plane, and a yellow density iso-contour of the droplets. 28 a)b)c) Figure 23: Schlieren images showing the evolution of the interaction of shock with R22 bubble in air at (a)t = 141.1 μs, (b) t = 240.5 μs and (c) t = 640.9 μs. 29 Density Pressure Velocity X Velocity Y Volume Fraction Schlieren Figure 24: Uniformly spaced time snapshots of Shock-induced R22 Bubble collapse in Water (SABW) at shock Mach number 1.70, with open boundary conditions on all sides (O), showing the evolution of pressure, density, x-velocity, y-velocity, volume fraction, and schlieren fromt 0 = 0.0 μs to t end = 7.0 μs. 30 A.4 Downsampling of Datasets Direct simulations on coarse grids are often of limited value because key phenomena, such as boundary layers or interface deformation, are poorly resolved, while training surrogate models on high-resolution data remains computationally challenging. A common workaround is to generate datasets at high resolution and then downsample them to a coarser grid by averaging primitive variables, such as velocity or pressure over neighboring cells. However, this naive averaging is not appropriate for compressible flows, since it does not preserve key conserved quantities such as momentum and energy. In this section, we present a conservative downsampling strategy that ensures conservation of total mass, momentum and energy when downsampling the dataset. A.4.1 Methodology When generating higher-resolution datasets with ALPACA, we also store the densities of each individual phase. Our methodology is designed to be general and applicable to any multiphase solver, which should, in principle, provide access to the densities of each phase. As all datasets in this work involve two phases, we save the phasic densities in addition to the mixed density field. The use of these phasic densities in the downsampling process is described in the following paragraphs. Another key property of multiphase simulation is the volume fraction,α, defining the share of each phase in a given cell. Here,α 1 defines the volume fraction of phase-1 andα 2 for phase-2 can be computed asα 2 = 1− α 1 . In cells where only one phase is present,αis either 0 or 1. However,α∈ (0, 1) in cells intersected by the interface. When downsampling the domain, these three equations need to be satisfied. The mass conservation states that the total mass in the domain must remain conserved during downsampling. Hence, the total mass in the fine grid cells must equal the total mass of the coarse grid cell, leading to Equation 6 for phase-1 (Φ 1 ), a similar equation is applicable for phase-2 (Φ 2 ) as well. δ 2 X i=1 ρ i,Φ 1 V i α i ! = ˆρ I,Φ 1 δ 2 X i=1 V i α i (6) The parameterδdefines the factor by which the grid is downsampled in each spatial direction and irepresents the index of the cells at the finer resolution. Equation 6 can be rearranged to compute the downsampled density of phase-1 in the coarse cellIas shown in Equation 7. Hereˆα I is the downsampled volume fraction of phase-1 in the coarse cellIand is computed as shown in Equation 8. ˆρ I,Φ 1 = P δ 2 i=1 ρ i,Φ 1 V i α i P δ 2 i=1 V i α i = 1 δ 2 ˆα I δ 2 X i=1 ρ i,Φ 1 α i (7) ˆα I = 1 δ 2 δ 2 X i=1 α i (8) To get the mixed density in the coarse cell, we combine the densities of both phases weighted by their volume fractions as shown in Equation 9. ˆρ I = ˆρ I,Φ 1 ˆα I + ˆρ I,Φ 2 (1− ˆα I )(9) Here we want to emphasize that both conservative and average downsampling lead to the same result for volume fraction and density fields. This is made evident by inserting Equation 7 into Equation 9. ˆρ I = ˆρ I,Φ 1 ˆα I + ˆρ I,Φ 2 (1− ˆα I ) = 1 δ 2 δ 2 X i=1 ρ i,Φ 1 α i + 1 δ 2 δ 2 X i=1 ρ i,Φ 2 (1− α i ) = 1 δ 2 δ 2 X i=1 ρ i (10) To get the downsampled x-velocity fieldu x , we use the momentum conservation in the x-direction. Similar to mass conservation, the total momentum in the fine grid cells must equal the total momentum 31 in the coarse grid cells as shown in Equation 11 for phase-1, and a similar equation applies for phase-2 and in every Cartesian direction. δ 2 X i=1 ρ i,Φ 1 u x,i,Φ 1 V i α i ! = ˆρ I,Φ 1 ˆu x,I,Φ 1 δ 2 X i=1 V i α i (11) In the above Equation,u x,i,Φ 1 is the x-component of the velocity of phase-1 for the fine grid celli. At the interface cells, we assume this quantity is equal for both phases and can therefore be written as u x,i . Outside the interface,u x,i =u x,i,Φ 1 α i +u x,i,Φ 2 (1− α i )=u x,i,Φ 1 forα i = 1. Thus Equation 11 can be rewritten as shown in Equation 12. δ 2 X i=1 ρ i,Φ 1 u x,i V i α i ! = ˆρ I,Φ 1 ˆu x,I,Φ 1 δ 2 X i=1 V i α i (12) From Equation 12, we get the downsampled x-velocity for phase-1 for the coarse grid cell I . ˆu x,I,Φ 1 = P δ 2 i=1 ρ i,Φ 1 u x,i V i α i ˆρ I,Φ 1 P δ 2 i=1 V i α i = 1 δ 2 ˆρ I,Φ 1 ˆα I δ 2 X i=1 ρ i,Φ 1 u x,i α i (13) To get the mixed momentum of the coarse cell, we combine the momentum of both phases weighted by their volume fractions as shown in Equation 14. Dividing by the mixed densityˆρ I gives us the downsampled mixed velocity in x-direction. ˆu x,I = ˆρ I,Φ 1 ˆu x,I,Φ 1 ˆα I + ˆρ I,Φ 2 ˆu x,I,Φ 2 (1− ˆα I ) ˆρ I (14) To compute the pressure fieldˆpand the total energy field ˆ Eof the downsampled flow, we need to use the stiffened gas EOS as shown in Equation 2, which is used by ALPACA for the generation of the datasets. By rearranging this equation we can compute the internal energy for phase-1 as shown in Equation 15. e i,Φ 1 = p i + γ Φ 1 p ∞,Φ 1 (γ Φ 1 − 1)ρ i,Φ 1 (15) Similar to the assumptions for no slip for the interfacial velocity, we assume that the interface is in mechanical equilibrium and thus pressure is equal for both phases at the interface cells such thatp i,Φ 1 =p i,Φ 2 =p i for these cells. Outside the interface,p i =p i,Φ 1 α i +p i,Φ 2 (1− α i )=p i,Φ 1 forα i = 1. Finally the total energy per unit mass for celliin phase-1 can be computed by combining the internal and kinetic energy as shown in Equation 16. E i,Φ 1 = e i,Φ 1 + 1 2 (u 2 x,i + u 2 y,i )(16) Here we use the mixed velocitiesu x,i andu y,i , because the velocities per phase are equal as described above. With Equation 16, we can now downsample the total energy field using the total energy conservation principle: The total energy in the fine grid cells must equal the total energy of the coarse grid cells as shown in Equation 17 for phase-1, and a similar equation applies for phase-2 as well. δ 2 X i=1 ρ i,Φ 1 E i,Φ 1 V i α i ! = ˆρ I,Φ 1 ˆ E I,Φ 1 δ 2 X i=1 V i α i (17) From this we get the downsampled total energy for phase-1 in coarse cellIas shown in Equation 18. ˆ E I,Φ 1 = P δ 2 i=1 ρ i,Φ 1 E i,Φ 1 V i α i ˆρ I,Φ 1 P δ 2 i=1 V i α i = 1 δ 2 ˆρ I,Φ 1 ˆα I δ 2 X i=1 ρ i,Φ 1 E i,Φ 1 α i (18) 32 Similarly to velocity, the mixed total energy can be computed as: ˆ E I = ˆρ I,Φ 1 ˆ E I,Φ 1 ˆα I + ˆρ I,Φ 2 ˆ E I,Φ 2 (1− ˆα I ) ˆρ I (19) We compute the internal energy per phase by rearranging Equation 16 and adapting it to the down- sampled values obtained from Equations 13 and 18. ˆe I,Φ 1 = ˆ E I,Φ 1 − 1 2 (ˆu 2 x,I,Φ 1 + ˆu 2 y,I,Φ 1 )(20) We refer to the EOS to compute the downsampled pressure for cellIin phase-1 as shown in Equation 21. ˆp I,Φ 1 = (γ Φ 1 − 1) ˆρ I,Φ 1 ˆe I,Φ 1 − γ Φ 1 p ∞,Φ 1 (21) To get the final downsampled pressureˆp I for the coarse cellI, we combine the pressures of both phases weighted by their volume fractions as shown in Equation 22. ˆp I = ˆp I,Φ 1 ˆα I + ˆp I,Φ 2 (1− ˆα I )(22) From these downsampled primitive fields we can further compute derived fields such as schlieren or vorticity. A.4.2 Results To compare the primitive and conservative downsampling methods, we compute the total mass, momentum and energy for the finer-resolution dataset at a particular timestep before downsampling in the entire domain as shown in Equation 23. Here N is the total number of cells in the domain. M tot = N X i=1 ρ i V i P tot,x/y/z = N X i=1 ρ i u x/y/z,i V i E tot = N X i=1 ρ i E i V i (23) After applying the average and conservative downsampling methods, we compute the total mass, momentum and energy again as shown in Equation 24. M tot,ds = M X I=1 ˆρ I ˆ V I P tot,ds,x/y/z = M X I=1 ˆρ I ˆu x/y/z,I ˆ V I E tot,ds = M X I=1 ˆρ I ˆ E I ˆ V I (24) In the above equation,M (= N/δ 2 )is the total number of coarse cells in the 2D downsampled domain. We choose a random trajectory from the 2D-shock-induced droplet breakup in air (SDBA) dataset and compare the results obtained from both downsampling methods. The average method downsamples each primitive field variable independently (e.g. density, velocity) by averaging the neighboring cells. The mass, momentum and energy of the averaged fields are computed after downsampling as shown in Figure 25. In contrast, the conservative method operates on the conserved quantities such as the momentum (ρuV) and energy (ρEV), ensuring conservation in the downsampled resolution and subsequently extracting the primitive variables of interest such as velocity. The absolute and relative errors in the conserved quantities obtained from conservative downsampling are in the order of machine precision. On the other hand, significant discrepancies are observed in the y-momentum, obtained from averaging the y-velocities on the fine grid - with relative errors ranging fromO(10 −2 ) -O(10), highlighting the need for a conservative downsampling approach. Although the relative error in the y-momentum during the initial timesteps is in the order ofO(10 4 ), this is primarily because the corresponding values on the original grid are nearly zero and division by such small numbers results in large relative errors. To gain insights into the impact of both downsampling methods, we visualize in Figure 26 the y-momentum field on the original resolution - (a) and the downsampled resolution - (b) and (c). When comparing the downsampled field to the original grid, the values obtained by averaging the y-velocities shows significant errors at the interface, whereas these errors are not present when using the conservative downsampling method. 33 0255075100125150 Time [ s] 0.111 0.112 0.113 0.114 0.115 0.116 0.117 Mass [kg] Total Mass Comparison Conservative Downsampling Average Downsampling Original Data 0255075100125150 Time [ s] 0 2 4 6 8 Momentum X [(kg m)/s] Total Momentum X Comparison Conservative Downsampling Average Downsampling Original Data 0255075100125150 Time [ s] 0.4 0.3 0.2 0.1 0.0 0.1 Momentum Y [(kg m)/s] Total Momentum Y Comparison Conservative Downsampling Average Downsampling Original Data 0255075100125150 Time [ s] 45000 46000 47000 48000 49000 50000 Total Energy [J] Total Energy Comparison Conservative Downsampling Average Downsampling Original Data 0255075100125150 Time [ s] 10 16 10 15 |Difference| [kg] Absolute Error in Mass Conservative Downsampling Average Downsampling 0255075100125150 Time [ s] 10 15 10 13 10 11 10 9 10 7 10 5 10 3 10 1 |Difference| [(kg m)/s] Absolute Error in Momentum X Conservative Downsampling Average Downsampling 0255075100125150 Time [ s] 10 15 10 13 10 11 10 9 10 7 10 5 10 3 10 1 |Difference| [(kg m)/s] Absolute Error in Momentum Y Conservative Downsampling Average Downsampling 0255075100125150 Time [ s] 10 14 10 11 10 8 10 5 10 2 10 1 |Difference| [J] Absolute Error in Total Energy Conservative Downsampling Average Downsampling 0255075100125150 Time [ s] 10 16 10 15 |Difference| / |Original| Relative Error in Mass Conservative Downsampling Average Downsampling 0255075100125150 Time [ s] 10 15 10 13 10 11 10 9 10 7 10 5 10 3 10 1 |Difference| / |Original| Relative Error in Momentum X Conservative Downsampling Average Downsampling 0255075100125150 Time [ s] 10 14 10 11 10 8 10 5 10 2 10 1 10 4 |Difference| / |Original| Relative Error in Momentum Y Conservative Downsampling Average Downsampling 0255075100125150 Time [ s] 10 15 10 13 10 11 10 9 10 7 10 5 10 3 |Difference| / |Original| Relative Error in Total Energy Conservative Downsampling Average Downsampling Shock Interaction with 5 Droplets in 2D Downsampling Factor: 4 - Original Resolution 512x512 Comparison of Downsampling Methods Figure 25: Comparison of mass, momentum and energy in the original grid [512x512] and downsam- pled grid [128x128] over time for the Shock-induced Droplet Breakup in Air (SDBA) case at shock Mach number 2.40 with 5 droplets. As observed in the y-momentum plot, downsampling by direct averaging leads to significant deviation from the reference, in contrast to the conservative method. 34 a)b)c)d) Figure 26: (a) y-momentum at the original resolution [512x512], (b) y-momentum at the downsampled resolution [128x128] obtained by averaging the y-velocities, and (c) y-momentum at the downsampled resolution [128x128] obtained by the conservative downsampling method. (d) Absolute difference between (b) and (c). The snapshot is taken at timet = 148 μs. Artifacts are observed at the interface when using the average downsampling method as seen in (b). 35 B Benchmarking Details B.1 Baseline Models We present additional details on the baseline models used in our benchmarking experiments. Our eval- uation covers a diverse set of state-of-the-art architectures trained from scratch with random initializa- tion, including convolutional models (CNeXt, CNO), spectral models (FFNO), and transformer-based architectures (ScOT). In addition, we fine-tune two pretrained models: Poseidon-B and DPOT-M. Details on the chosen hyperparameters for the baselines are as follows: 1.Convolutional Neural Operator (CNO): A CNO, much like a U-Net, processes an input function through a sequence of encoder layers that progressively decrease spatial resolution while increasing channel depth. This encoded representation is then passed to decoder layers, which reverse the process by restoring spatial resolution and reducing the number of channels. At corresponding spatial scales (or spectral bands), encoder and decoder features are connected via ResNet-style skip connections. The central idea of the architecture is its upsampling and downsampling scheme, where low-pass filtering is applied to avoid introducing artificial high-frequency components during resampling. This approach aligns with the Shannon sampling theorem, ensuring that discrete representations remain consistent with the underlying continuous functions. Table 5: CNO hyperparameters. Hyperparameter50M Latent channels256 Depth24 # blocks (except at neck)8 # blocks (at neck)14 Channel multiplier24 Activationcustom LeakyReLU 2.ConvNeXt (CNeXt): A ConvNeXt is a modern convolutional neural network architecture introduced to bridge the performance gap between traditional CNNs and vision transformers. It builds on the classical ResNet design but incorporates several transformer-inspired updates, such as large-kernel depthwise convolutions, Layer Normalization, GELU activations, and a patch-based stem. Organized as a hierarchical model with multiple stages that progressively reduce spatial resolution while increasing channel capacity, ConvNeXt serves as a powerful and efficient feature extractor for a wide range of vision tasks. Despite remaining fully convolutional, its design achieves performance comparable to transformer-based models. The 3D datasets are trained on the 10M variant and the 2D datasets are trained on the 50M variant. Table 6: CNeXt hyperparameters. Hyperparameter10M50M Latent channels3280 # stages44 # blocks per stage (except at neck)11 # blocks (at neck)11 Channel multiplier22 ActivationGELUGELU 3.Factorized Fourrier Neural Operator(FFNO): Like FNO, it operates in the Fourier domain to capture global interactions, but instead of learning a full multi-dimensional spectral kernel, FFNO factorizes the operator into a sequence of one-dimensional Fourier transformations applied along each spatial dimension. This separable formulation signifi- cantly reduces the number of parameters and computational cost while maintaining strong 36 expressive power. Similar to ConvNeXt, the 3D datasets are trained using the 10M variant, while the 2D datasets are trained using the 50M variant. Table 7: FFNO hyperparameters. Hyperparameter10M50M Latent channels96240 FNO layers812 Modes2016 MLP expansion factor44 # MLP layers22 Activation in MLPReLUReLU Activation in input and output headsGELUGELU 4.Scalable Operator Transformer (ScOT) and Poseidon: The ScOT model is built upon the Poseidon framework (Herde et al., 2024), where Poseidon provides the foundation- model framework for learning generalizable solution operators, and ScOT serves as its core architectural backbone. At its core, ScOT employs a hierarchical transformer design inspired by vision transformers, using a window-based mechanism in which the input domain is divided into a uniform grid of non-overlapping patches, with added support for non-square inputs. Each patch is averaged via a shared spatial weighting scheme and then linearly projected into a latent embedding space defined by the number of latent channels, yielding a piecewise-constant latent representation that reduces the cost of global attention while retaining essential local structure. This embedded representation is then processed through multiple stages of SwinV2 Transformer blocks (Liu et al., 2021) arranged in a UNet-like hierarchy, where feature maps are progressively downsampled and later upsampled, with the number of blocks per stage governed by the “depths” hyperparameter. Within each stage, attention is computed using windowed multi-head self-attention (MHSA), restricting interactions to local regions for efficiency, while shifted windows between layers enable cross-window communication and capture global context. The encoder incorporates patch merging to decrease spatial resolution and increase feature dimensionality for learning global patterns, whereas the decoder restores resolution via patch expansion, with skip connections implemented through ConvNeXt blocks (Liu et al., 2022)—controlled by the “skip-connections” hyperparameter—linking corresponding encoder and decoder stages to preserve fine-grained information. Table 8: ScOT hyperparameters. Hyperparameter50M Latent channels150 Patch size4 Depths[4, 4, 4] # attention heads[6, 12, 24] Skip connections[3, 3, 0] Window size16 MLP ratio4.0 ActivationGELU Table 9: Poseidon-B hyperparameters. Hyperparameter158M Latent channels96 Patch size4 Depths[8, 8, 8, 8] # attention heads[3, 6, 12, 24] Skip connections[2, 2, 2, 0] Window size16 MLP ratio4.0 ActivationGELU 5.Denoising Pre-trained Operator Transformer (DPOT): DPOT is a large-scale neural operator architecture designed to learn solution mappings of PDEs. It is built around a Fourier transformer backbone, where attention is performed in the frequency domain to efficiently capture global spatial dependencies, enabling the model to approximate integral operators over function spaces. The architecture combines temporal aggregation layers to encode information from multiple time steps with multi-head Fourier attention layers, which apply learnable transformations in Fourier space to model complex multi-scale dynamics. Furthermore, DPOT is paired with an auto-regressive denoising pre-training strategy, where the model predicts future states from noise-corrupted inputs to improve robustness and generalization. This design allows DPOT to scale to very large models and 37 datasets, functioning as a foundation model for PDEs that can be fine-tuned across diverse downstream tasks. Table 10: DPOT-M hyperparameters. Hyperparameter122M Latent channels1024 Patch size8 # blocks8 MLP ratio4.0 Depth12 Modes32 Mixing TypeAFNO Time AggregationExponential MLP ActivationGELU B.2 Training Protocol As mentioned in the main text, we employ a many-to-one training strategy for the models intialized with random weights, where 4 historic snapshots are passed to the model as input and the next snapshot is predicted and compared with the target. We use this setup for benchmarking the baselines mentioned in the previous section, but the benchmarking repository is capabale of handling many- to-many scenarios as well. For the pre-trained models, we adhere to the original input–output configurations specified in their respective works. A key central feature of our benchmarking repository is the ability to integrate conditioning parameters into all baseline architectures via custom normalization layers (Figure 27) where intermediate feature statistics are modulated through learned scale and shift terms derived from the conditioning vector. Consider an input to the custom norm,x∈R B×S×L , whereBdenotes the batch size,Srepresent the spatial resolution, andLis the latent dimension. The custom norm is similar to a standard LayerNorm (Ba et al., 2016), with the modification that the affine parameters (γ θ andβ θ ) are made functions of the conditioning parameter vectorc(Perez et al., 2017), wherec∈R B×F , andFis the number of conditioning features. Figure 27: Illustration of the AdaNorm conditioning module, applicable to all baselines.Fdenotes the number of conditioning parameters, and L the latent dimension. 38 y(c) = γ θ (c)· x− μ σ + β θ (c) μ = 1 BSL B X b=1 S X s=1 L X ℓ=1 x σ = v u u t 1 BSL B X b=1 S X s=1 L X ℓ=1 (x− μ) 2 + ε. (25) In our benchmarking repository, we have several other conditioning methods implemented, but in this work we focus on Adaptive Normalization (AdaNorm), which is a simple and lightweight mapping of the conditioning feature vector c, through an MLP to the desired latent dimension L. The hyperparameters used for training the 2D and 3D baselines are shown in Table 11. Table 11: Training hyperparameters used across experiments. Hyperparameter2D3D Train time time-restricted to 24 hours time-restricted to 24 hours on a 4 NVIDIA-A100 on a single NVIDIA-H100 Checkpointing metric RMSE RMSE for baselines Batch size161 Validation split ratio0.20.2 Optimizer AdamW AdamW Learning rate (LR) 1e-3 1e-4 LR scheduler cosine cosine Weight decay 1e-6 1e-6 B.3 Evaluation Protocol Model performance is evaluated on a held-out test set. During inference, each model is applied autoregressively, starting from the initial condition and rolling out for the full length of a test trajectory. Inference produces prediction and target tensors of shape(N,R,T,C,S), whereNis the number of test trajectories,Rthe number of rollout steps, andTthe number of predicted timesteps per rollout. Each metric is reduced over the spatial dimensions only, yielding an error tensor of shape E ∈R N×(RT )×C . This preserves the dependence on trajectory, prediction time, and channel. The resulting tensor enables evaluation at multiple levels. First, the mean and standard deviation over the trajectory dimension are computed asμ,σ ∈R (RT )×C , which describe the temporal evolution of the error and the variability across test trajectories, reported separately for each channel. A scalar overall loss for each metricμ overall across all trajectories of the test dataset is obtained by aggregating over the remaining dimensions ofμ, i.e. over prediction time and channels, in accordance with the formulation used by that metric. Tables 22 through 51 showcase this scalar overall loss for all the metrics in this work for different datasets. B.4 Metrics As mentioned in the main text, we perform two sets of ablations from a metric perspective. In the first set of experiments, we perform training of the baselines with only MSE loss and in the second set of experiments, we train the baselines with a composite loss which is a combination of several metrics: MSE, H1, SSIM and IRMSE. The initial weights for these metrics is specified in Table 12. During inference, we evaluate the performance of the trained baselines on a broader set of metrics, as mentioned in Table 13. 39 Table 12: Initial loss coefficients for the composite loss. Initial loss weights default λ MSE 1.0 λ H1 2.0 λ SSIM 0.5 λ IRMSE 1.0 Table 13: Overview of implemented evaluation metrics. Metrics are grouped by conceptual type. NameAcronymTypeSource Mean squared errorMSEPointwise Mean absolute errorMAEPointwise Variance-normalized RMSEVRMSEPointwise Ohanaetal. (2024) H 1 -seminormH1Pointwise Hartmann et al. (2019); Cho et al. (2024a) H 2 -seminormH2PointwiseHartmann et al. (2019); Cho et al. (2024a) Interface-focused RMSEIRMSEFeature-focusedHassanetal. (2023, 2025) Shock-focused RMSESRMSEFeature-focusedCooper- Baldock et al. (2024) Multilevel waveletMLWSpectralPrantletal. (2022) Wavelet frequency-binned RMSEwfRMSESpectralTakamoto et al. (2022);Prantl et al. (2022) Structural similarity indexSSIMStructure-awareNilssonand Akenine-Möller (2020) Pearson correlation coefficientPCCStructure-awareHu et al. (2026); Maoetal. (2025) Integral quantities of interest:iqRMSEPhysics-informedTakamoto et al. (2022) Mass, Momentum, Kinetic energy, Enstrophy, Center of mass displacement The inputs to each loss function are the prediction (y) and target (ˆy) tensors, where ˆy,y ∈R B×T×C×S B, T, C, S correspond to the batch, time, channel and spatial dimensions respectively. Spatial sums are written over the full discrete grid, and metrics are computed per sample before averaging over the batch dimension. For readability, the batch index is omitted in the formulas described in this section. B.4.1 Pointwise metrics The Mean Square Error (MSE) and its varient the Root Mean Square Error (RMSE) are defined as shown in Equations 26 and 27 respectively. L MSE = 1 N T,C,S X T,C,S (ˆy− y) 2 (26) 40 L RMSE = s 1 N T,C,S X T,C,S (ˆy− y) 2 (27) The Varience scaled Root Mean Square Error (VRMSE) metric divides the mean squared error by the variance of the target field as shown in Equation 28. L VRMSE = v u u t 1 N T,C X T,C P S (ˆy− y) 2 P S (y− ̄y) 2 + ε (28) where ̄y is the spatial mean of the target field and ε is 10 −8 . Another useful variant of the MSE family is the Normalized Root Mean Square Error (nRMSE), which normalizes the RMSE by the mean squared magnitude of the target field. The nRMSE instead normalizes by the mean squared magnitude of the reference field as shown in Equation 29. L nRMSE = v u u t 1 N T,C X T,C P S (ˆy− y) 2 P S y 2 ,(29) The Sobolev loss of order K is defined as L Sob,K = λ 0 1 N T,C,S X T,C,S (ˆy− y) 2 + K X k=1 λ k 1 N T,C,S X T,C,S ∇ (k) S ˆy−∇ (k) S y 2 ,(30) where∇ (k) denotes thek-th order spatial derivative,λ k are weighting coefficients, and∥·∥ 2 is the pointwise squared Euclidean norm of the derivative vector/tensor Hartmann et al. (2019); Cho et al. (2024a). In this work, the derivative terms of the Sobolev loss are each defined as separate metrics, corresponding to the H 1 and H 2 components of the Sobolev space. TheH 1 semi-norm metric, as shown in Figure 28 penalizes mismatches in first-order spatial deriva- tives as defined in Equation 31: L H 1 = 1 N T,C,S X T,C,S d X j=1 ∂ x j ˆy− ∂ x j y 2 .(31) The H 2 semi-norm penalizes discrepancies in second-order derivatives as shown in Equation 32. L H 2 = 1 N T,C,S X T,C,S d X i=1 d X j=1 ∂ 2 x i x j ˆy− ∂ 2 x i x j y 2 .(32) Spatial derivatives are computed on regular grids using fixed convolution stencils, supporting both standard central-difference and Sobel kernels, the latter combining differentiation with mild local smoothing. B.4.2 Feature-focused metrics The Interface-focused Root Mean Square Error (IRMSE) restricts error evaluation to the vicinity of the phase boundary (Figure 29) as defined in Equation 33. L IRMSE = v u u t 1 N TC X TC P S m (ˆy− y) 2 P S m + ε .(33) Here m denotes the mask of spatial locations corresponding to the interface region. In the original formulation,mis defined as the zero level-set of a signed-distance fieldφ(x). In this work, the interface is identified using a threshold on the density fieldρ. Since a hard threshold would 41 Figure 28: H 1 semi-norm metric. introduce discontinuities in the loss, a soft interface mask is constructed using sigmoid transitions around the density bounds ρ min and ρ max , which maintains differentiability: m = σ ρ− ρ min s σ ρ max − ρ s ,(34) where s controls the softness of the transition. Figure 29: IRMSE metric. Table 14: Configuration settings for interface-focused RMSE. Hyperparamter SABW SDBA SRBA Threshold field density density density Threshold range (350, 550) (350, 550) (4.5, 8.5) Blur sigma 3.0 3.0 3.0 Threshold softness 0.05 0.05 0.05 The GRMSE emphasizes regions where the target field exhibits strong spatial variation as described in Cooper-Baldock et al. (2024). We define a member of this family called shock-RMSE with a special focus in capturing strong shocks from the pressure field. A spatial weighting functionw(y)is constructed from the magnitude of the gradient of the target pressure field: w ∝∥∇ S y∥ 2 ,(35) 42 The Shock-focused Root Mean Square Error (SRMSE) is then defined as shown in Equation 36. L SRMSE = v u u t 1 N T,C X T,C P S w (ˆy− y) 2 P S w + ε .(36) In this formulation, spatial locations with large gradients receive higher weights. For this imple- mentation, we definewusing the spatial gradients of the pressure field (Figure 30). This effectively produces a mask which emphasizes shock fronts, providing a complementary metric to IRMSE. Figure 30: SRMSE metric. Table 15: Configuration settings for shock-focused RMSE. Hyperparameter default Gradient threshold field pressure Gradient threshold range (0.3, 4.0) Threshold softness 0.03 Blur sigma 1.5 Derivative stencil Central Difference Both IRMSE and SRMSE optionally include Gaussian smoothing, as proposed in the original GRMSE formulation Cooper-Baldock et al. (2024). A Gaussian filter with standard deviationσ is applied to the weighting masks before computing the loss. This is implemented as separable convolution, where a one-dimensional Gaussian kernel is applied sequentially along each spatial axis. The smoothing reduces high-frequency noise in the weighting fields and leads to more stable gradients during training. B.4.3 Spectral metrics In this work, spectral decomposition is implemented using the multi-level DWT. The key idea is to represent a field in terms of its spatially localized frequency content across different spatial scales. At each decomposition leveli, the input field is split into a low-frequency component and a high-frequency component, LF i (y), HF i (y),(37) and analogously for the prediction ˆy, LF i (ˆy), HF i (ˆy).(38) The low-frequency componentLF i contains the coarse, slowly varying structure of the field, while the high-frequency componentHF i captures finer-scale detail. The transform is then applied recursively 43 only to the low-frequency component. Repeating this procedure overIlevels yields one final coarse approximation LF I together with a hierarchy of detail coefficients, HF 1 ,HF 2 ,...,HF I ,(39) whereHF 1 corresponds to the finest details and higher levels represent progressively coarser spatial scales (Figure 31). Figure 31: Conceptual illustration of the DWT, recursively decomposing an inputyinto frequency- band contributions. For multidimensional fields, the transform is applied independently along each spatial dimension. This produces directional detail subbands (e.g., horizontal, vertical, and diagonal features in 2D). Collectively, these directional components form the set of high-frequency coefficientsHF i at leveli. In this work, orthogonal Daubechies wavelets are used, following their usage in related studies Prantl et al. (2022). For the Multilevel Wavelet (MLW) metric, the DWT is first computed over the spatial dimensions of yandˆy, yielding the detail coefficientsHF i at each levelias depicted in Figure 32. The coefficients are then passed through a logarithmic transform: w i (y) = log 2 |HF i (y)| + ε , w i (ˆy) = log 2 |HF i (ˆy)| + ε ,(40) where ε is a small constant for numerical stability. The wavelet spectral loss is then defined as: L w = X i ∥w i (y)− w i (ˆy)∥ 1 ,(41) with the sum taken over all decomposition levels and orientations. The logarithmic scaling makes the more loss sensitive to discrepancies in fine-scale structure (i.e., thin interfaces and small-scale oscillations). When applied to trajectories with multiple timesteps, the loss optionally includes a separate 1D wavelet transform over the time dimension, applied independently at each spatial location. In this case, the loss is constructed as the sum of spatial and temporal contributions, where α and β are the respective weighting coefficients: L MLW = αL ws + βL wt (42) The settings used for α and β are provided in Table 16. Table 16: Configuration settings for MultilevelWaveletLoss. HyperparameterDefault Wavelet db2 Alpha 1 Beta 0.1 DWT mode reflect The Wavelet frequency-binned Root Mean Square Error (wfRMSE) metric measures RMSE across wavelet frequency bands. This metric is inspired by the frequency-binned Fourier RMSE used in Takamoto et al. (2022), but replaces global Fourier modes with localized wavelet bands to avoid periodicity assumptions. 44 Figure 32: Multilevel wavelet loss (MLW) metric. For each wavelet level i, the wfRMSE of the detail coefficients is computed as RMSE i = q HF i (ˆy)− HF i (y) 2 .(43) The total wavelet frequency-binned RMSE is defined as a weighted sum over scales, L wfRMSE = X i α i RMSE i ,(44) where the weightsα i control the relative contribution of different frequency bands. By adjustingα i , the metric can be tuned to prioritize specific ranges of dynamically important scales. Table 17: Configuration settings for wfRMSE. Hyperparameter high mid low Wavelet db2 db2 db2 Levels 3 3 3 Level weights (1,0,0) (0,1,0) (0,0,1) DWT mode reflect reflect reflect B.4.4 Structure-aware metrics The Pearson Correlation Coefficient (PCC) metric measures the linear correlation between two signals and is defined as: r = 1 N T,C X T,C P S (y− ̄y)(ˆy− ̄ ˆy) q P S (y− ̄y) 2 q P S (ˆy− ̄ ˆy) 2 .(45) Here ̄ ˆy and ̄y denote the means of the predictions and labels. The corresponding loss is defined as: L PCC = 1− r.(46) PCC is invariant under affine transformations of the inputs, making it insensitive to uniform shifts or global scaling differences between fields. Minimizing1− rtherefore encourages the prediction to reproduce the overall covariance structure of the target, promoting spatial alignment of features. The Structural Similarity Index (SSIM) for two given fieldsy, ˆyis evaluated over sliding spatial windows and decomposes similarity into three components: luminance, contrast, and structure as described in Nilsson and Akenine-Möller (2020) and Venkataramanan et al. (2021), SSIM(y, ˆy) = l(y, ˆy)c(y, ˆy)s(y, ˆy),(47) 45 These components are defined as: l(y, ˆy) = 2μ 1 μ 2 + C 1 μ 2 1 + μ 2 2 + C 1 , c(y, ˆy) = 2σ 1 σ 2 + C 2 σ 2 1 + σ 2 2 + C 2 , s(y, ˆy) = σ 12 + C 3 σ 1 σ 2 + C 3 ,(48) whereμ i ,σ 2 i , andσ 12 denote the local means, variances, and covariance of the two fields computed within a spatial window.C 1 ,C 2 , andC 3 are small constants that stabilize division (typically C 3 = C 2 /2). The global SSIM score is obtained by spatially averaging local scores: MSSIM(y, ˆy) = 1 N W X w∈W SSIM(y, ˆy;w).(49) The resulting values range from−1to1, with1denoting perfect similarity. When used as a training objective, this metric is written in the form (Figure 33): L SSIM = 1− MSSIM.(50) Figure 33: SSIM metric. The hyperparameters used for the SSIM metric are provided in Table 18. Table 18: Configuration settings for SSIM. Hyperparameter default Window size 11 Sigma 1.5 K1 0.01 K2 0.03 L 10 B.4.5 Physics-informed metrics Letc(y)denote a scalar quantity of interest. On a discrete grid, we approximate the integral of this quantity over Ω by a sum over grid cells. Q(y) = X S c(y) ∆V(51) where ∆V is the cell volume or quadrature weight. The Domain-integrated Quantity of Interest RMSE (iqRMSE) is defined as: 46 Table 19: Quantities included in the iqRMSE metric. QuantityExpression Mass R Ω ρdΩ Momentum along axis k R Ω ρu k dΩ Kinetic energy R Ω 1 2 ρ|u| 2 dΩ Enstrophy R Ω 1 2 |∇× u| 2 dΩ Center of mass along axis k R Ω ρx k dΩ R Ω ρdΩ iqRMSE = v u u t 1 N T X T Q(ˆy)− Q(y) Q(y) + ε 2 .(52) This metric measures discrepancies in the total amount of the conserved quantity, independent of how errors are distributed spatially within the domain (Figure 34). The implemented quantities are summarized in Table 19. In the implementation, all predicted and target fields are first converted back to physical units before computing the metric. To keep the resulting values well-scaled, the conservation error is normalized by the corresponding target quantity, resulting in an nRMSE-style normalization that yields loss values of order O(1). Figure 34: iqRMSE metric. B.5 Adaptive loss balancing framework Given a set of M loss components: L i M i=1 , the scalar training objective is written as J (θ) = M X i=1 λ i L i (θ), λ i > 0,(53) where λ i denotes the weight assigned to loss component i. A generalized loss-balancing framework is illustrated in Figure 35. The scalarized training lossJ (θ) is computed via Equation 53, and the constituent loss components are individually sampled and stored. The same is optionally performed for component gradients, whereby only aggregated statistics across parameters are stored (indicated here byg j ) as opposed to full gradient tensors, in order to 47 reduce memory overhead. These histories form the input to the loss weighting strategy. The sampling frequency and update frequency are configurable independently. By default, loss values or gradients are sampled at every training step, while the loss weights are updated once per epoch. This makes it possible to trade off responsiveness against computational overhead. In this work, we employ two loss balancing strategies: SoftAdapt as shown in algorithm 1 and the Gradient Normalization (GradNorm) as shown in algorithm 2. Figure 35: Conceptual illustration of the loss balancing framework. B.5.1 SoftAdapt SoftAdapt Heydari et al. (2019) adjusts the weights of loss components according to their recent rate of change. The key idea is to prioritize objectives whose loss has decreased more slowly (or increased), indicating that the model is currently struggling to optimize them. To estimate this behavior, SoftAdapt approximates the recent slope of each loss component using a short history of past values. LetL (t) i denote the value of loss componentiat training stept. The recent rate of change is approximated ass (t) i = L (t) i −L (t−1) i . These slopes are converted into loss component weights through a softmax transformation λ (t) i = exp(βs (t) i ) P T j=1 exp(βs (t) j ) ,(54) whereβis a tunable temperature parameter that controls how strongly the weighting emphasizes poorly improving losses. Forβ > 0, losses with larger positive slopes receive higher weights, while β < 0favors losses that are decreasing fastest. Settingβ = 0results in uniform weighting. The hyperparameter value chosen for this work is presented in Table 20. Table 20: Configuration settings for SoftAdapt. Hyperparameter default Temperature (β) 1.0 B.5.2 Gradient Normalization (GradNorm) GradNorm (Chen et al., 2017) updates the loss component weights such that the gradient norm of each weighted loss component matches a target determined by the component’s relative training speed. In practice, the method is applied to a shared parameter subsetθ sh ⊂ θ, typically the last shared layer. For taski, the gradient norm is defined asG (i) θ = ∥∇ θ (λ i L i )∥ 2 . The target gradient norm ˆ G (i) θ sh = G θ sh [r i ] α is based on the normalized loss ratio ̃ L i and the relative inverse training rate r i , r i = ̃ L i 1 M P j ̃ L j , ̃ L i = L i L i (0) ,(55) 48 Algorithm 1 Training with SoftAdapt 1: Choose temperature parameter β 2: Initialize weights λ i ← 1 M ∀i 3: Store initial lossesL (0) i 4: for each epoch do 5:for each training step do 6:Train network with weighted lossJ = P M i=1 λ i L i 7:Record current lossesL (t) i 8:end for 9:Compute slope estimates s (t) i =L (t) i −L (t−1) i 10:Compute weights λ (t) i (Equation 54) 11:Normalize weights such that P M i=1 λ i = 1 12: end for and the asymmetry hyperparameterα. The hyperparameter values selected for this work are presented in Table 21. Algorithm 2 Training with GradNorm 1: Initialize task weights λ i ← λ i (0)∀i 2: Store initial task lossesL i (0) 3: for each epoch do 4:for each training step do 5:Train network for with weighted lossJ = P M i=1 λ i L i 6:Compute gradient norms on shared parameters θ sh as G (i) θ sh =∥∇ θ sh (λ i L i )∥ 2 7:end for 8:Compute average gradient norm across the M loss components: G θ sh = 1 M P M i=1 G (i) θ sh 9:Compute normalized losses ̃ L i and inverse training rates r i (Equation 55) 10:Define target gradient norms ˆ G (i) θ sh = G θ sh [r i ] α 11:Compute GradNorm objectiveL grad = P M i=1 G (i) θ sh − ˆ G (i) θ sh 12:Update task weights using∇ λ i L grad (targets ˆ G (i) θ sh treated as constants) 13:Renormalize weights such that P M i=1 λ i = M 14: end for Table 21: Configuration settings for GradNorm. Hyperparameter default Asymmetry parameter (α) 1.0 Loss weight learning rate 0.001 B.6 Computational Resources Each trajectory of the high fidelity 2D-datasets and the 3D-SDBA (SSOOSS) datasets was generated on an HPC-cluster using a single Intel-Xeon(R) Platinum 8480+ node with 112 CPU-cores, while a single 3D-SABW (SSOOSS) trajectory was generated on a single node of Intel Skylake Xeon Platinum 8174 with 24 CPU-cores. For training surrogate models, we employ a Data Distributed Parallel strategy on 4 NVIDIA-A100 GPUs each with 40GB of VRAM for the 2D baselines, while the 3D baselines were trained on 4 NVIDIA-H100 GPUs, each with 90GB VRAM. 49 C Detailed results C.1 Metric performance comparison Figure 36 shows the difference in the prediction quality for the last timestep of a 2D-SDBA trajectory using FFNO with two different training strategies and also comparison with the ground truth. It can be observed that the predictions having better SSIM values have more structural coherence with the ground truth even though RMSE indicates equally good predictions for both. Ground Truth a) b) c) Prediction (GradNorm-Training) Prediction (MSE-Training) Density Density Density Figure 36: (a) Reference density snapshot for a randomly chosen 2D-SDBA dataset at t=70 s. (b) Corresponding FFNO (50M) prediction trained with composite loss and further using GradNorm adaptive loss weighting strategy. RMSE = 0.46 and SSIM = 0.28.(c) Corresponding FFNO prediction trained with only MSE loss. RMSE = 0.46 and SSIM = 0.32 As the results in the main text indicate clear advantages of using composite losses with adaptive loss-weighting strategies, a natural question concerns the additional computational cost incurred by these methods. We examine this trade-off from two perspectives: the peak GPU memory requirement, shown in Figure 37, and the training time required to complete one epoch, shown in Figure 38. Since GradNorm requires gradient information to adapt the loss weights, it is expected to introduce additional memory and computational overhead. In contrast, SoftAdapt does not require gradient storage and therefore introduces negligible overhead in both memory usage and training time per epoch. C.2 Inference metric evolution over rollout steps This section presents the evolution of the inference metrics across rollout steps for the test dataset trajectories. The metrics are averaged across all trajectories in the test dataset at each rollout step and the corresponding mean and standard deviation are presented in the following figures. The figures are organized by dataset, with each figure corresponding to a specific dataset and containing the results for all baselines and training strategies evaluated on that dataset. The metrics are categorized into pointwise, spectral, feature-focused, structure-aware, and physics-informed. Poseidon-M has a different starting point in all the plots shown in this section, since the first prediction timestep depends on the number of input timesteps required by the baseline and the stride between those timesteps. 50 CNO-50M FFNO-50M ConvNeXt-50M ScOT-50M DPOT-M-122M Poseidon-B-158M 0.0% 2.0% 4.0% 6.0% 8.0% 10.0% Percentage increase relative to MSE training loss 5.1% 2.9% 10.0% 5.0% 1.1% 5.8% 0.0%0.0%0.0%0.0% 1.1% 0.1% Percentage Increase in Maximum GPU Memory Allocation Averaged across 2D Datasets for Composite Training Loss Relative to MSE Training Loss GradNorm SoftAdapt Figure 37: Peak GPU memory allocation during training CNO-50M FFNO-50M ConvNeXt-50M ScOT-50M DPOT-M-122M Poseidon-B-158M 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 Epoch duration [min] Epoch Duration Averaged across 2D Datasets for Different Training Strategies MSE SoftAdapt GradNorm Figure 38: Epoch duration 51 C.2.1 2D shock-induced air bubble collapse in water with open boundaries [2D-SABW (O)] 255075100 Timestep 10 3 10 2 10 1 H1 Training with MSE Loss 255075100 Timestep 10 3 10 2 10 1 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 4 10 3 10 2 10 1 Training with Composite Loss (GradNorm) 255075100 Timestep 10 3 10 2 H2 255075100 Timestep 10 3 10 2 255075100 Timestep 10 3 10 2 255075100 Timestep 10 3 10 2 10 1 10 0 MAE 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 2 10 1 10 0 RMSE 255075100 Timestep 10 2 10 1 10 0 255075100 Timestep 10 2 10 1 10 0 255075100 Timestep 10 1 10 0 10 1 10 2 10 3 VRMSE 255075100 Timestep 10 2 10 1 10 0 10 1 10 2 255075100 Timestep 10 2 10 1 10 0 10 1 10 2 10 3 Evolution of Pointwise Rollout Metrics for 2D Shock Induced Air Bubble Collapse in Water [2D-SABW (O)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 39: Pointwise evaluation of baselines on the 2D-SABW (O) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields for each timestep. 52 255075100 Timestep 10 3 10 2 10 1 wfRMSE-high Training with MSE Loss 255075100 Timestep 10 2 10 1 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 2 10 1 Training with Composite Loss (GradNorm) 255075100 Timestep 10 2 10 1 wfRMSE-mid 255075100 Timestep 10 3 10 2 10 1 255075100 Timestep 10 3 10 2 10 1 255075100 Timestep 10 2 10 1 10 0 wfRMSE-low 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 0 3 × 10 1 4 × 10 1 6 × 10 1 2 × 10 0 MLW 255075100 Timestep 10 0 3 × 10 1 4 × 10 1 6 × 10 1 255075100 Timestep 10 0 3 × 10 1 4 × 10 1 6 × 10 1 Evolution of Spectral Rollout Metrics for 2D Shock Induced Air Bubble Collapse in Water [2D-SABW (O)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 40: Wavelet based spectral metrics for the 2D-SABW (O) test dataset, averaged across all fields for each timestep. 53 20406080100 Timestep 10 2 10 1 10 0 IRMSE Training with MSE Loss 20406080100 Timestep 10 2 10 1 10 0 Training with Composite Loss (SoftAdapt) 20406080100 Timestep 10 2 10 1 10 0 Training with Composite Loss (GradNorm) 20406080100 Timestep 10 2 10 1 10 0 SRMSE 20406080100 Timestep 10 3 10 2 10 1 10 0 20406080100 Timestep 10 2 10 1 10 0 Evolution of Feature-Focused Rollout Metrics for 2D Shock Induced Air Bubble Collapse in Water [2D-SABW (O)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 41: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 2D-SABW (O) test dataset, averaged across all fields for each timestep. 20406080100 Timestep 10 3 10 2 10 1 10 0 PCC Training with MSE Loss 20406080100 Timestep 10 4 10 3 10 2 10 1 10 0 Training with Composite Loss (SoftAdapt) 20406080100 Timestep 10 4 10 3 10 2 10 1 10 0 Training with Composite Loss (GradNorm) 20406080100 Timestep 10 1 10 0 SSIM 20406080100 Timestep 10 2 10 1 10 0 20406080100 Timestep 10 2 10 1 10 0 Evolution of Structure-Aware Rollout Metrics for 2D Shock Induced Air Bubble Collapse in Water [2D-SABW (O)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 42: Structure-aware metrics for the 2D-SABW (O) test dataset, averaged across all fields for each timestep. 54 255075100 Timestep 10 3 10 2 10 1 10 0 iqRMSE-mass Training with MSE Loss 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 3 10 2 10 1 10 0 Training with Composite Loss (GradNorm) 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 iqRMSE-momentum_x 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 2 10 1 10 0 10 1 255075100 Timestep 10 0 10 1 10 2 10 3 iqRMSE-momentum_y 255075100 Timestep 10 0 10 1 10 2 10 3 255075100 Timestep 10 0 10 1 10 2 10 3 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 10 1 iqRMSE-kinetic-energy 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 3 10 2 10 1 10 0 iqRMSE-enstrophy 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 Evolution of Physics-Informed Rollout Metrics for 2D Shock Induced Air Bubble Collapse in Water [2D-SABW (O)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 43: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 2D-SABW (O) test dataset, averaged across all fields for each timestep. 55 C.2.2 2D shock-induced air bubble collapse in water with symmetry boundaries [2D-SABW (SSOO)] 255075100 Timestep 10 5 10 4 10 3 10 2 10 1 H1 Training with MSE Loss 255075100 Timestep 10 4 10 3 10 2 10 1 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 3 10 2 10 1 Training with Composite Loss (GradNorm) 255075100 Timestep 10 3 10 2 H2 255075100 Timestep 10 4 10 3 10 2 255075100 Timestep 10 4 10 3 10 2 255075100 Timestep 10 3 10 2 10 1 10 0 MAE 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 2 10 1 10 0 RMSE 255075100 Timestep 10 2 10 1 10 0 255075100 Timestep 10 2 10 1 10 0 255075100 Timestep 10 2 10 1 10 0 10 1 10 2 10 3 VRMSE 255075100 Timestep 10 1 10 0 10 1 10 2 10 3 255075100 Timestep 10 1 10 0 10 1 10 2 10 3 Evolution of Pointwise Rollout Metrics for 2D Shock Induced Air Bubble Collapse in Water [2D-SABW (SSOO)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 44: Pointwise evaluation of baselines on the 2D-SABW (SSOO) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields for each timestep. 56 255075100 Timestep 10 3 10 2 10 1 wfRMSE-high Training with MSE Loss 255075100 Timestep 10 2 10 1 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 2 10 1 Training with Composite Loss (GradNorm) 255075100 Timestep 10 2 10 1 wfRMSE-mid 255075100 Timestep 10 3 10 2 10 1 255075100 Timestep 10 2 10 1 255075100 Timestep 10 2 10 1 10 0 wfRMSE-low 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 2 10 1 10 0 255075100 Timestep 10 0 3 × 10 1 4 × 10 1 6 × 10 1 2 × 10 0 MLW 255075100 Timestep 10 0 3 × 10 1 4 × 10 1 6 × 10 1 255075100 Timestep 10 0 3 × 10 1 4 × 10 1 6 × 10 1 Evolution of Spectral Rollout Metrics for 2D Shock Induced Air Bubble Collapse in Water [2D-SABW (SSOO)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 45: Wavelet based spectral metrics for the 2D-SABW (SSOO) test dataset, averaged across all fields for each timestep. 57 20406080100 Timestep 10 1 10 0 IRMSE Training with MSE Loss 20406080100 Timestep 10 2 10 1 10 0 Training with Composite Loss (SoftAdapt) 20406080100 Timestep 10 1 10 0 Training with Composite Loss (GradNorm) 20406080100 Timestep 10 2 10 1 10 0 SRMSE 20406080100 Timestep 10 2 10 1 10 0 20406080100 Timestep 10 2 10 1 10 0 Evolution of Feature-Focused Rollout Metrics for 2D Shock Induced Air Bubble Collapse in Water [2D-SABW (SSOO)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 46: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 2D-SABW (SSOO) test dataset, averaged across all fields for each timestep. 20406080100 Timestep 10 4 10 3 10 2 10 1 10 0 PCC Training with MSE Loss 20406080100 Timestep 10 5 10 4 10 3 10 2 10 1 10 0 Training with Composite Loss (SoftAdapt) 20406080100 Timestep 10 3 10 2 10 1 10 0 Training with Composite Loss (GradNorm) 20406080100 Timestep 10 1 10 0 SSIM 20406080100 Timestep 10 2 10 1 10 0 20406080100 Timestep 10 1 10 0 Evolution of Structure-Aware Rollout Metrics for 2D Shock Induced Air Bubble Collapse in Water [2D-SABW (SSOO)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 47: Structure-aware metrics for the 2D-SABW (SSOO) test dataset, averaged across all fields for each timestep. 58 255075100 Timestep 10 5 10 4 10 3 10 2 10 1 10 0 10 1 iqRMSE-mass Training with MSE Loss 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 3 10 2 10 1 10 0 Training with Composite Loss (GradNorm) 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 10 1 10 2 iqRMSE-momentum_x 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 1 10 2 10 3 iqRMSE-momentum_y 255075100 Timestep 10 1 10 2 10 3 255075100 Timestep 10 1 10 2 10 3 255075100 Timestep 10 3 10 1 10 1 10 3 iqRMSE-kinetic-energy 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 3 10 2 10 1 10 0 iqRMSE-enstrophy 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 Evolution of Physics-Informed Rollout Metrics for 2D Shock Induced Air Bubble Collapse in Water [2D-SABW (SSOO)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 48: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 2D-SABW (SSOO) test dataset, averaged across all fields for each timestep. 59 C.2.3 2D shock-induced R22 bubble collapse in air with open boundaries [2D-SRBA (O)] 255075100 Timestep 10 4 10 3 10 2 10 1 H1 Training with MSE Loss 255075100 Timestep 10 3 10 2 10 1 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 3 10 2 10 1 Training with Composite Loss (GradNorm) 255075100 Timestep 10 3 10 2 H2 255075100 Timestep 10 3 10 2 255075100 Timestep 10 6 10 5 10 4 10 3 10 2 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 MAE 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 255075100 Timestep 10 2 10 1 10 0 RMSE 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 2 10 1 10 0 255075100 Timestep 10 2 10 1 10 0 10 1 10 2 VRMSE 255075100 Timestep 10 2 10 1 10 0 10 1 10 2 255075100 Timestep 10 2 10 1 10 0 10 1 10 2 Evolution of Pointwise Rollout Metrics for 2D Shock Induced R22 Bubble Collapse in Air [2D-SRBA (O)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 49: Pointwise evaluation of baselines on the 2D-SRBA (O) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields for each timestep. 60 255075100 Timestep 10 3 10 2 10 1 wfRMSE-high Training with MSE Loss 255075100 Timestep 10 2 10 1 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 2 10 1 Training with Composite Loss (GradNorm) 255075100 Timestep 10 2 10 1 wfRMSE-mid 255075100 Timestep 10 3 10 2 10 1 255075100 Timestep 10 3 10 2 10 1 255075100 Timestep 10 2 10 1 10 0 wfRMSE-low 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 0 3 × 10 1 4 × 10 1 6 × 10 1 MLW 255075100 Timestep 10 0 3 × 10 1 4 × 10 1 6 × 10 1 255075100 Timestep 10 0 3 × 10 1 4 × 10 1 6 × 10 1 Evolution of Spectral Rollout Metrics for 2D Shock Induced R22 Bubble Collapse in Air [2D-SRBA (O)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 50: Wavelet based spectral metrics for the 2D-SRBA (O) test dataset, averaged across all fields for each timestep. 61 20406080100 Timestep 10 4 10 3 10 2 10 1 10 0 IRMSE Training with MSE Loss 20406080100 Timestep 10 4 10 3 10 2 10 1 10 0 Training with Composite Loss (SoftAdapt) 20406080100 Timestep 10 2 10 1 10 0 Training with Composite Loss (GradNorm) 20406080100 Timestep 10 2 10 1 10 0 SRMSE 20406080100 Timestep 10 3 10 2 10 1 10 0 20406080100 Timestep 10 3 10 2 10 1 10 0 Evolution of Feature-Focused Rollout Metrics for 2D Shock Induced R22 Bubble Collapse in Air [2D-SRBA (O)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 51: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 2D-SRBA (O) test dataset, averaged across all fields for each timestep. 20406080100 Timestep 10 4 10 3 10 2 10 1 10 0 PCC Training with MSE Loss 20406080100 Timestep 10 4 10 3 10 2 10 1 10 0 Training with Composite Loss (SoftAdapt) 20406080100 Timestep 10 4 10 3 10 2 10 1 10 0 Training with Composite Loss (GradNorm) 20406080100 Timestep 10 1 10 0 SSIM 20406080100 Timestep 10 2 10 1 10 0 20406080100 Timestep 10 2 10 1 10 0 Evolution of Structure-Aware Rollout Metrics for 2D Shock Induced R22 Bubble Collapse in Air [2D-SRBA (O)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 52: Structure-aware metrics for the 2D-SRBA (O) test dataset, averaged across all fields for each timestep. 62 255075100 Timestep 10 3 10 2 10 1 10 0 iqRMSE-mass Training with MSE Loss 255075100 Timestep 10 2 10 1 10 0 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 3 10 2 10 1 10 0 Training with Composite Loss (GradNorm) 255075100 Timestep 10 2 10 1 10 0 10 1 iqRMSE-momentum_x 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 2 10 4 10 6 10 8 iqRMSE-momentum_y 255075100 Timestep 10 2 10 4 10 6 10 8 255075100 Timestep 10 2 10 4 10 6 10 8 255075100 Timestep 10 2 10 1 10 0 10 1 10 2 iqRMSE-kinetic-energy 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 10 2 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 2 10 1 10 0 iqRMSE-enstrophy 255075100 Timestep 10 3 10 2 10 1 10 0 255075100 Timestep 10 4 10 3 10 2 10 1 10 0 Evolution of Physics-Informed Rollout Metrics for 2D Shock Induced R22 Bubble Collapse in Air [2D-SRBA (O)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 53: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 2D-SRBA (O) test dataset, averaged across all fields for each timestep. 63 C.2.4 2D shock-induced droplet breakup in air with symmetry boundaries [2D-SDBA (SSOO)] 255075100 Timestep 10 3 10 2 10 1 H1 Training with MSE Loss 255075100 Timestep 10 3 10 2 10 1 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 3 10 2 10 1 Training with Composite Loss (GradNorm) 255075100 Timestep 10 3 10 2 H2 255075100 Timestep 10 3 10 2 255075100 Timestep 10 3 10 2 255075100 Timestep 10 1 10 0 MAE 255075100 Timestep 10 1 10 0 255075100 Timestep 10 1 10 0 255075100 Timestep 10 1 10 0 RMSE 255075100 Timestep 10 1 10 0 255075100 Timestep 10 1 10 0 255075100 Timestep 10 1 10 0 10 1 VRMSE 255075100 Timestep 10 1 10 0 10 1 255075100 Timestep 10 1 10 0 10 1 Evolution of Pointwise Rollout Metrics for 2D Shock Induced Droplet Breakup in Air [2D-SDBA (SSOO)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 54: Pointwise evaluation of baselines on the 2D-SDBA (SSOO) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields for each timestep. 64 255075100 Timestep 10 1 3 × 10 2 4 × 10 2 6 × 10 2 wfRMSE-high Training with MSE Loss 255075100 Timestep 10 1 2 × 10 2 3 × 10 2 4 × 10 2 6 × 10 2 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 1 2 × 10 2 3 × 10 2 4 × 10 2 6 × 10 2 Training with Composite Loss (GradNorm) 255075100 Timestep 10 1 6 × 10 2 2 × 10 1 3 × 10 1 4 × 10 1 wfRMSE-mid 255075100 Timestep 10 1 255075100 Timestep 10 1 255075100 Timestep 10 1 10 0 wfRMSE-low 255075100 Timestep 10 1 10 0 255075100 Timestep 10 1 10 0 255075100 Timestep 10 0 3 × 10 1 4 × 10 1 6 × 10 1 MLW 255075100 Timestep 10 0 2 × 10 1 3 × 10 1 4 × 10 1 6 × 10 1 255075100 Timestep 10 0 2 × 10 1 3 × 10 1 4 × 10 1 6 × 10 1 Evolution of Spectral Rollout Metrics for 2D Shock Induced Droplet Breakup in Air [2D-SDBA (SSOO)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 55: Wavelet based spectral metrics for the 2D-SDBA (SSOO) test dataset, averaged across all fields for each timestep. 65 20406080100 Timestep 10 0 IRMSE Training with MSE Loss 20406080100 Timestep 10 0 Training with Composite Loss (SoftAdapt) 20406080100 Timestep 10 0 Training with Composite Loss (GradNorm) 20406080100 Timestep 10 1 10 0 SRMSE 20406080100 Timestep 10 1 10 0 20406080100 Timestep 10 1 10 0 Evolution of Feature-Focused Rollout Metrics for 2D Shock Induced Droplet Breakup in Air [2D-SDBA (SSOO)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 56: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 2D-SDBA (SSOO) test dataset, averaged across all fields for each timestep. 20406080100 Timestep 10 3 10 2 10 1 PCC Training with MSE Loss 20406080100 Timestep 10 2 10 1 Training with Composite Loss (SoftAdapt) 20406080100 Timestep 10 3 10 2 10 1 Training with Composite Loss (GradNorm) 20406080100 Timestep 10 1 SSIM 20406080100 Timestep 10 1 20406080100 Timestep 10 1 Evolution of Structure-Aware Rollout Metrics for 2D Shock Induced Droplet Breakup in Air [2D-SDBA (SSOO)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 57: Structure-aware metrics for the 2D-SDBA (SSOO) test dataset, averaged across all fields for each timestep. 66 255075100 Timestep 10 2 10 1 10 0 10 1 iqRMSE-mass Training with MSE Loss 255075100 Timestep 10 2 10 1 10 0 10 1 Training with Composite Loss (SoftAdapt) 255075100 Timestep 10 2 10 1 10 0 10 1 Training with Composite Loss (GradNorm) 255075100 Timestep 10 2 10 1 10 0 10 1 iqRMSE-momentum_x 255075100 Timestep 10 2 10 1 10 0 10 1 255075100 Timestep 10 2 10 1 10 0 10 1 255075100 Timestep 10 0 10 1 10 2 iqRMSE-momentum_y 255075100 Timestep 10 0 10 1 10 2 255075100 Timestep 10 0 10 1 10 2 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 iqRMSE-kinetic-energy 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 3 10 2 10 1 10 0 10 1 255075100 Timestep 10 1 10 0 10 1 iqRMSE-enstrophy 255075100 Timestep 10 1 10 0 10 1 255075100 Timestep 10 1 10 0 10 1 255075100 Timestep 10 3 10 2 10 1 iqRMSE-com_x 255075100 Timestep 10 4 10 3 10 2 10 1 255075100 Timestep 10 3 10 2 10 1 255075100 Timestep 10 4 10 3 10 2 10 1 iqRMSE-com_y 255075100 Timestep 10 4 10 3 10 2 10 1 255075100 Timestep 10 5 10 4 10 3 10 2 10 1 Evolution of Physics-Informed Rollout Metrics for 2D Shock Induced Droplet Breakup in Air [2D-SDBA (SSOO)] CNO (50M) ConvNeXt (50M) DPOT-M (122M) FFNO (50M) Poseidon-B (158M) ScOT (50M) Figure 58: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 2D-SDBA (SSOO) test dataset, averaged across all fields for each timestep. 67 C.2.5 3D shock-induced air bubble collapse in water with symmetry boundaries [3D-SABW (SSOOSS)] 20304050 Timestep 10 2 10 1 H1 Training with MSE Loss 20304050 Timestep 10 1 2 × 10 2 3 × 10 2 4 × 10 2 6 × 10 2 Training with Composite Loss (SoftAdapt) 20304050 Timestep 10 1 2 × 10 2 3 × 10 2 4 × 10 2 6 × 10 2 Training with Composite Loss (GradNorm) 20304050 Timestep 10 2 10 1 H2 20304050 Timestep 10 1 20304050 Timestep 10 2 10 1 20304050 Timestep 10 2 10 1 MAE 20304050 Timestep 10 2 10 1 20304050 Timestep 10 2 10 1 20304050 Timestep 10 1 RMSE 20304050 Timestep 10 1 20304050 Timestep 10 1 20304050 Timestep 10 1 10 0 10 1 10 2 VRMSE 20304050 Timestep 10 0 10 1 20304050 Timestep 10 0 10 1 Evolution of Pointwise Rollout Metrics for 3D Shock Induced Air Bubble Collapse in Water [3D-SABW (SSOOSS)] ConvNeXt (10M) FFNO (10M) Figure 59: Pointwise evaluation of baselines on the 3D-SABW (SSOOSS) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields for each timestep. 68 20304050 Timestep 10 1 2 × 10 2 3 × 10 2 4 × 10 2 6 × 10 2 wfRMSE-high Training with MSE Loss 20304050 Timestep 10 1 2 × 10 2 3 × 10 2 4 × 10 2 6 × 10 2 Training with Composite Loss (SoftAdapt) 20304050 Timestep 10 1 2 × 10 2 3 × 10 2 4 × 10 2 6 × 10 2 Training with Composite Loss (GradNorm) 20304050 Timestep 10 1 wfRMSE-mid 20304050 Timestep 10 1 20304050 Timestep 10 1 20304050 Timestep 10 1 10 0 wfRMSE-low 20304050 Timestep 10 1 10 0 20304050 Timestep 10 1 10 0 20304050 Timestep 10 0 6 × 10 1 2 × 10 0 3 × 10 0 4 × 10 0 MLW 20304050 Timestep 10 0 6 × 10 1 2 × 10 0 3 × 10 0 4 × 10 0 20304050 Timestep 10 0 6 × 10 1 2 × 10 0 3 × 10 0 4 × 10 0 Evolution of Spectral Rollout Metrics for 3D Shock Induced Air Bubble Collapse in Water [3D-SABW (SSOOSS)] ConvNeXt (10M) FFNO (10M) Figure 60: Wavelet based spectral metrics for the 3D-SABW (SSOOSS) test dataset, averaged across all fields for each timestep. 69 20304050 Timestep 10 0 IRMSE Training with MSE Loss 20304050 Timestep 10 1 10 0 Training with Composite Loss (SoftAdapt) 20304050 Timestep 10 1 10 0 Training with Composite Loss (GradNorm) 20304050 Timestep 10 1 10 0 SRMSE 20304050 Timestep 10 1 10 0 20304050 Timestep 10 1 10 0 Evolution of Feature-Focused Rollout Metrics for 3D Shock Induced Air Bubble Collapse in Water [3D-SABW (SSOOSS)] ConvNeXt (10M) FFNO (10M) Figure 61: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 3D-SABW (SSOOSS) test dataset, averaged across all fields for each timestep. 20304050 Timestep 10 2 10 1 PCC Training with MSE Loss 20304050 Timestep 10 2 10 1 Training with Composite Loss (SoftAdapt) 20304050 Timestep 10 2 10 1 Training with Composite Loss (GradNorm) 20304050 Timestep 2 × 10 1 3 × 10 1 4 × 10 1 6 × 10 1 SSIM 20304050 Timestep 10 1 20304050 Timestep 10 1 2 × 10 1 3 × 10 1 4 × 10 1 6 × 10 1 Evolution of Structure-Aware Rollout Metrics for 3D Shock Induced Air Bubble Collapse in Water [3D-SABW (SSOOSS)] ConvNeXt (10M) FFNO (10M) Figure 62: Structure-aware metrics for the 3D-SABW (SSOOSS) test dataset, averaged across all fields for each timestep. 70 20304050 Timestep 10 3 10 2 10 1 10 0 10 1 10 2 iqRMSE-mass Training with MSE Loss 20304050 Timestep 10 4 10 3 10 2 Training with Composite Loss (SoftAdapt) 20304050 Timestep 10 3 10 2 10 1 10 0 10 1 10 2 Training with Composite Loss (GradNorm) 20304050 Timestep 10 3 10 2 10 1 10 0 10 1 10 2 iqRMSE-momentum_x 20304050 Timestep 10 2 10 1 20304050 Timestep 10 3 10 2 10 1 10 0 10 1 10 2 20304050 Timestep 10 1 10 3 10 5 10 7 iqRMSE-momentum_y 20304050 Timestep 10 1 10 3 10 5 10 7 20304050 Timestep 10 0 10 2 10 4 10 6 10 8 20304050 Timestep 10 3 10 2 10 1 10 0 10 1 10 2 iqRMSE-kinetic-energy 20304050 Timestep 10 2 10 1 20304050 Timestep 10 3 10 2 10 1 10 0 10 1 10 2 20304050 Timestep 10 3 10 2 10 1 10 0 iqRMSE-enstrophy 20304050 Timestep 10 2 10 1 10 0 20304050 Timestep 10 2 10 1 10 0 Evolution of Physics-Informed Rollout Metrics for 3D Shock Induced Air Bubble Collapse in Water [3D-SABW (SSOOSS)] ConvNeXt (10M) FFNO (10M) Figure 63: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 3D-SABW (SSOOSS) test dataset, averaged across all fields for each timestep. 71 C.2.6 3D shock-induced droplet-breakup in air with symmetry boundaries [3D-SDBA (SSOOSS)] 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 H1 Training with MSE Loss 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 Training with Composite Loss (SoftAdapt) 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 Training with Composite Loss (GradNorm) 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 H2 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 MAE 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 RMSE 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 VRMSE 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 Evolution of Pointwise Rollout Metrics for 3D Shock Induced Droplet Breakup in Air [3D-SDBA (SSOOSS)] ConvNeXt (10M) FFNO (10M) Figure 64: Pointwise evaluation of baselines on the 3D-SDBA (SSOOSS) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields for each timestep. 72 406080100 Timestep 10 1 10 0 wfRMSE-high Training with MSE Loss 406080100 Timestep 10 1 10 0 10 1 Training with Composite Loss (SoftAdapt) 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 Training with Composite Loss (GradNorm) 406080100 Timestep 10 1 10 0 wfRMSE-mid 406080100 Timestep 10 1 10 0 10 1 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 406080100 Timestep 10 0 10 1 wfRMSE-low 406080100 Timestep 10 1 10 0 10 1 10 2 406080100 Timestep 10 2 10 1 10 0 10 1 10 2 406080100 Timestep 10 0 2 × 10 0 MLW 406080100 Timestep 10 0 406080100 Timestep 10 0 Evolution of Spectral Rollout Metrics for 3D Shock Induced Droplet Breakup in Air [3D-SDBA (SSOOSS)] ConvNeXt (10M) FFNO (10M) Figure 65: Wavelet based spectral metrics for the 3D-SDBA (SSOOSS) test dataset, averaged across all fields for each timestep. 73 406080100 Timestep 10 0 10 1 IRMSE Training with MSE Loss 406080100 Timestep 10 0 10 1 10 2 Training with Composite Loss (SoftAdapt) 406080100 Timestep 10 1 10 0 10 1 10 2 Training with Composite Loss (GradNorm) 406080100 Timestep 10 1 10 0 10 1 SRMSE 406080100 Timestep 10 1 10 0 10 1 406080100 Timestep 10 1 10 0 10 1 10 2 Evolution of Feature-Focused Rollout Metrics for 3D Shock Induced Droplet Breakup in Air [3D-SDBA (SSOOSS)] ConvNeXt (10M) FFNO (10M) Figure 66: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 3D-SDBA (SSOOSS) test dataset, averaged across all fields for each timestep. 406080100 Timestep 10 2 10 1 PCC Training with MSE Loss 406080100 Timestep 10 2 10 1 Training with Composite Loss (SoftAdapt) 406080100 Timestep 10 1 Training with Composite Loss (GradNorm) 406080100 Timestep 3 × 10 1 4 × 10 1 6 × 10 1 SSIM 406080100 Timestep 3 × 10 1 4 × 10 1 5 × 10 1 406080100 Timestep 3 × 10 1 4 × 10 1 6 × 10 1 Evolution of Structure-Aware Rollout Metrics for 3D Shock Induced Droplet Breakup in Air [3D-SDBA (SSOOSS)] ConvNeXt (10M) FFNO (10M) Figure 67: Structure-aware metrics for the 3D-SDBA (SSOOSS) test dataset, averaged across all fields for each timestep. 74 406080100 Timestep 10 0 10 1 10 2 iqRMSE-mass Training with MSE Loss 406080100 Timestep 10 0 10 1 10 2 Training with Composite Loss (SoftAdapt) 406080100 Timestep 10 0 10 1 10 2 Training with Composite Loss (GradNorm) 406080100 Timestep 10 0 10 1 10 2 10 3 iqRMSE-momentum_x 406080100 Timestep 10 0 10 1 10 2 10 3 406080100 Timestep 10 0 10 1 10 2 10 3 406080100 Timestep 10 0 10 1 10 2 10 3 iqRMSE-momentum_y 406080100 Timestep 10 0 10 1 10 2 10 3 406080100 Timestep 10 0 10 1 10 2 10 3 406080100 Timestep 10 0 10 1 10 2 iqRMSE-kinetic-energy 406080100 Timestep 10 0 10 1 10 2 406080100 Timestep 10 0 10 1 10 2 406080100 Timestep 10 2 10 1 10 0 iqRMSE-enstrophy 406080100 Timestep 10 1 10 0 406080100 Timestep 10 1 10 0 406080100 Timestep 10 1 2 × 10 1 3 × 10 1 4 × 10 1 iqRMSE-com_x 406080100 Timestep 10 1 4 × 10 2 6 × 10 2 406080100 Timestep 10 2 10 1 406080100 Timestep 10 2 10 1 iqRMSE-com_y 406080100 Timestep 10 4 10 3 10 2 10 1 406080100 Timestep 2 × 10 2 3 × 10 2 4 × 10 2 6 × 10 2 Evolution of Physics-Informed Rollout Metrics for 3D Shock Induced Droplet Breakup in Air [3D-SDBA (SSOOSS)] ConvNeXt (10M) FFNO (10M) Figure 68: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 3D-SDBA (SSOOSS) test dataset, averaged across all fields for each timestep. 75 C.3 Overall inference metrics The following tables report the long-term rollout inference metrics. For each trajectory in the test dataset, the trained baseline model is applied autoregressively, starting from the initial condition and continuing until the end of the trajectory. The evolution of the metrics, averaged across all the trajectories, for each rollout for the test dataset is presented in section C.2. The final metric presented in this section is obtained by averaging the metric values across all trajectories and rollout steps. The tables are organized by dataset, with each table corresponding to a specific dataset and containing the results for all baselines and training strategies evaluated on that dataset. The metrics are categorized into pointwise, spectral, feature-focused, structure-aware, and physics-informed. Table 22: Pointwise evaluation of baselines on the 2D-SABW (O) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields and all rollout steps. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossH1H2MAERMSEVRMSE ConvNeXt (50M) MSE loss 0.02520.00680.06400.22550.2858 Comp.loss (SA) 0.02510.00670.06910.22990.2841 Comp.loss (GN) 0.02440.00700.07190.22540.2875 CNO (50M) MSE loss 0.02700.00730.09510.25750.3655 Comp.loss (SA)0.02770.00690.09510.26180.3447 Comp.loss (GN) 0.02370.00640.07770.22650.2967 FFNO (50M) MSE loss0.02300.00630.05710.20810.2389 Comp.loss (SA) 0.02270.00620.05800.20520.2320 Comp.loss (GN)0.02310.00620.06410.21290.2466 ScOT (50M) MSE loss 0.02300.00720.05560.20300.2589 Comp.loss (SA) 0.02700.00770.07730.24880.3022 Comp.loss (GN)0.02460.00740.07170.22320.2918 DPOT-M (122M) MSE loss0.01760.00570.04210.16350.1956 Comp.loss (SA)0.01980.00610.05350.18500.2305 Comp.loss (GN)0.02010.00630.05690.18600.2336 Poseidon-B (158M) MSE loss0.05160.01000.50880.8666>10 Comp.loss (SA) 0.05010.01000.43360.7727>10 Comp.loss (GN)0.05140.01000.47410.8279>10 Table 23: Wavelet based spectral metrics for the 2D-SABW (O) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining losswfRMSE-wfRMSE-wfRMSE-MLW highmidlow ConvNeXt (50M) MSE loss0.06570.18130.43260.8191 Comp.loss (SA)0.06440.17840.43060.7089 Comp.loss (GN)0.06890.17640.41920.7489 CNO (50M) MSE loss0.06950.19340.45780.6822 Comp.loss (SA) 0.06770.19470.47220.5831 Comp.loss (GN)0.06470.17950.42330.5763 FFNO (50M) MSE loss 0.05890.16090.39190.7298 Comp.loss (SA)0.05840.15960.38560.5833 Comp.loss (GN)0.05980.16370.39810.5956 ScOT (50M) MSE loss0.06780.16790.37980.7929 Comp.loss (SA)0.07230.18710.43310.6881 Comp.loss (GN)0.07340.18320.40800.7610 DPOT-M (122M) MSE loss 0.05790.14440.31350.7733 Comp.loss (SA)0.06170.15950.35130.6953 Comp.loss (GN)0.06330.16200.35290.7456 Poseidon-B (158M) MSE loss0.08550.28710.85371.1695 Comp.loss (SA)0.08510.27770.82221.0185 Comp.loss (GN)0.08530.28310.84541.0685 76 Table 24: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 2D-SABW (O) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossIRMSESRMSE ConvNeXt (50M) MSE loss1.06770.3405 Comp.loss (SA)1.04350.3442 Comp.loss (GN) 1.02790.3385 CNO (50M) MSE loss1.08740.3775 Comp.loss (SA)1.10150.3827 Comp.loss (GN)0.98460.3392 FFNO (50M) MSE loss0.97910.3104 Comp.loss (SA)0.94940.3068 Comp.loss (GN) 0.96650.3164 ScOT (50M) MSE loss0.98360.3082 Comp.loss (SA)1.05470.3559 Comp.loss (GN) 0.98410.3357 DPOT-M (122M) MSE loss 0.80760.2496 Comp.loss (SA)0.91010.2855 Comp.loss (GN)0.86970.2891 Poseidon-B (158M) MSE loss 1.77781.0315 Comp.loss (SA)1.71290.9423 Comp.loss (GN) 1.75270.9931 Table 25: Structure-aware metrics for the 2D-SABW (O) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossPCCSSIM ConvNeXt (50M) MSE loss0.05310.1938 Comp.loss (SA) 0.05600.1560 Comp.loss (GN)0.05430.1736 CNO (50M) MSE loss0.07820.2512 Comp.loss (SA)0.07630.2042 Comp.loss (GN)0.05680.1841 FFNO (50M) MSE loss0.04930.1505 Comp.loss (SA)0.04800.1375 Comp.loss (GN)0.05050.1537 ScOT (50M) MSE loss 0.04440.1754 Comp.loss (SA)0.06150.1702 Comp.loss (GN)0.05360.1846 DPOT-M (122M) MSE loss0.02770.1327 Comp.loss (SA)0.03720.1328 Comp.loss (GN)0.03750.1535 Poseidon-B (158M) MSE loss0.47570.4742 Comp.loss (SA) 0.44450.4145 Comp.loss (GN)0.47000.4379 77 Table 26: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 2D-SABW (O) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossiqRMSE-iqRMSE-iqRMSE-iqRMSE-iqRMSE- massmomentummomentumkinetic-energyenstrophy _x_y ConvNeXt (50M) MSE loss0.00840.0122>100.01490.1232 Comp.loss (SA)0.01620.0269>100.03170.1205 Comp.loss (GN)0.02550.0350>100.04080.1198 CNO (50M) MSE loss0.03550.0753>100.07130.1142 Comp.loss (SA) 0.03030.0491>100.06010.1271 Comp.loss (GN)0.01880.0284>100.02980.1046 FFNO (50M) MSE loss0.00410.0060>100.00930.0829 Comp.loss (SA)0.00710.0112>100.01180.0752 Comp.loss (GN)0.01450.0208>100.02240.0717 ScOT (50M) MSE loss 0.00870.0132>100.02660.0996 Comp.loss (SA)0.01260.0163>100.02620.1457 Comp.loss (GN) 0.01930.0237>100.03530.1647 DPOT-M (122M) MSE loss0.00800.0167>100.02150.0525 Comp.loss (SA)0.00940.0174>100.02660.0736 Comp.loss (GN) 0.01350.0218>100.02370.0672 Poseidon-B (158M) MSE loss 0.45282.6884>106.16791.2095 Comp.loss (SA)0.33941.7529>103.46450.9530 Comp.loss (GN) 0.38402.0922>104.37531.0224 Table 27: Pointwise evaluation of baselines on the 2D-SABW (SSOO) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields and all rollout steps. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossH1H2MAERMSEVRMSE ConvNeXt (50M) MSE loss0.02650.00740.07260.22950.4207 Comp.loss (SA) 0.02600.00700.07530.22890.3504 Comp.loss (GN)0.02480.00720.07890.22500.3923 CNO (50M) MSE loss 0.02810.00740.09990.25170.5123 Comp.loss (SA)0.02840.00710.09980.25930.3719 Comp.loss (GN)0.02450.00660.08650.22860.3434 FFNO (50M) MSE loss0.02400.00670.06400.20890.2620 Comp.loss (SA)0.02360.00660.06650.20750.2307 Comp.loss (GN)0.02440.00660.07410.21820.2843 ScOT (50M) MSE loss0.02440.00790.05860.20350.3706 Comp.loss (SA)0.02520.00750.06570.21580.3335 Comp.loss (GN)0.02520.00790.07760.22410.3944 DPOT-M (122M) MSE loss0.01930.00630.05670.17610.2650 Comp.loss (SA)0.02050.00640.06220.18750.2645 Comp.loss (GN)0.02080.00660.06500.18960.2754 Poseidon-B (158M) MSE loss0.05410.01070.47300.8278>10 Comp.loss (SA) 0.05530.01110.51470.8799>10 Comp.loss (GN)0.05390.01070.51760.8806>10 78 Table 28: Wavelet based spectral metrics for the 2D-SABW (SSOO) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining losswfRMSE-wfRMSE-wfRMSE-MLW highmidlow ConvNeXt (50M) MSE loss0.06940.19040.44520.8141 Comp.loss (SA)0.06670.18420.43880.6925 Comp.loss (GN)0.07050.18110.42110.7396 CNO (50M) MSE loss0.07090.19650.46420.6835 Comp.loss (SA) 0.06910.19520.47660.5666 Comp.loss (GN)0.06700.18410.43210.5659 FFNO (50M) MSE loss0.06190.16810.40120.7349 Comp.loss (SA)0.06090.16610.39460.5846 Comp.loss (GN)0.06350.17350.41310.6009 ScOT (50M) MSE loss0.07120.17370.38300.7634 Comp.loss (SA) 0.07160.17980.40480.6736 Comp.loss (GN) 0.07630.18770.40860.7326 DPOT-M (122M) MSE loss0.06330.15870.33650.7750 Comp.loss (SA)0.06350.16220.35090.6991 Comp.loss (GN) 0.06470.16470.35790.7301 Poseidon-B (158M) MSE loss 0.08980.29720.87161.1604 Comp.loss (SA)0.09150.29610.88961.0478 Comp.loss (GN) 0.08910.29180.88121.0689 Table 29: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 2D-SABW (SSOO) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossIRMSESRMSE ConvNeXt (50M) MSE loss 1.15290.3310 Comp.loss (SA)1.11360.3287 Comp.loss (GN) 1.06450.3216 CNO (50M) MSE loss 1.13790.3523 Comp.loss (SA)1.15550.3631 Comp.loss (GN)1.01930.3266 FFNO (50M) MSE loss 1.04270.2973 Comp.loss (SA)1.00460.2935 Comp.loss (GN)1.01380.3083 ScOT (50M) MSE loss1.05160.2933 Comp.loss (SA)1.07200.3138 Comp.loss (GN)1.04380.3203 DPOT-M (122M) MSE loss0.90330.2568 Comp.loss (SA)0.92640.2714 Comp.loss (GN)0.92460.2749 Poseidon-B (158M) MSE loss 1.89230.9896 Comp.loss (SA) 1.97111.0434 Comp.loss (GN)1.98291.0453 79 Table 30: Structure-aware metrics for the 2D-SABW (SSOO) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossPCCSSIM ConvNeXt (50M) MSE loss0.04720.2037 Comp.loss (SA)0.04900.1630 Comp.loss (GN)0.04740.1822 CNO (50M) MSE loss 0.06120.2649 Comp.loss (SA)0.06630.2041 Comp.loss (GN)0.05270.1967 FFNO (50M) MSE loss0.04280.1596 Comp.loss (SA)0.04240.1499 Comp.loss (GN) 0.04550.1709 ScOT (50M) MSE loss0.03820.1666 Comp.loss (SA)0.04390.1487 Comp.loss (GN)0.04700.1870 DPOT-M (122M) MSE loss0.02900.1669 Comp.loss (SA)0.03280.1448 Comp.loss (GN)0.03470.1637 Poseidon-B (158M) MSE loss0.45210.4427 Comp.loss (SA) 0.48680.4518 Comp.loss (GN)0.47770.4537 Table 31: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 2D-SABW (SSOO) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossiqRMSE-iqRMSE-iqRMSE-iqRMSE-iqRMSE- massmomentummomentumkinetic-energyenstrophy _x_y ConvNeXt (50M) MSE loss 0.00820.0133>100.01800.1171 Comp.loss (SA)0.01720.0263>100.03170.1267 Comp.loss (GN) 0.01940.0277>100.03280.1330 CNO (50M) MSE loss 0.62552.7342>107.48990.1267 Comp.loss (SA)0.01960.0405>100.04580.1075 Comp.loss (GN)0.01410.0214>100.02440.1015 FFNO (50M) MSE loss0.00370.0058>100.00930.0842 Comp.loss (SA)0.00650.0089>100.01290.0783 Comp.loss (GN)0.01140.0147>100.01640.0711 ScOT (50M) MSE loss1.02672.6316>106.00670.0976 Comp.loss (SA)0.01100.0137>100.01850.1271 Comp.loss (GN)0.02110.0248>100.03570.1740 DPOT-M (122M) MSE loss 0.02320.0503>100.06510.0566 Comp.loss (SA)0.00970.0237>100.03650.0811 Comp.loss (GN)0.01940.0203>100.02110.0720 Poseidon-B (158M) MSE loss0.39132.3253>105.16861.1436 Comp.loss (SA)0.43022.6383>106.13921.2832 Comp.loss (GN)0.43982.7058>106.40821.2329 80 Table 32: Pointwise evaluation of baselines on the 2D-SRBA (O) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields and all rollout steps. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossH1H2MAERMSEVRMSE ConvNeXt (50M) MSE loss0.02940.01110.05090.18470.1979 Comp.loss (SA) 0.02830.01070.04930.17770.1833 Comp.loss (GN)0.02770.01060.05100.17730.1849 CNO (50M) MSE loss0.03650.01240.07560.23290.2538 Comp.loss (SA) 0.03140.01110.06350.20080.2091 Comp.loss (GN)0.03380.01160.06960.21480.2293 FFNO (50M) MSE loss0.02700.01040.04430.16740.1667 Comp.loss (SA)0.02570.01000.04140.15860.1557 Comp.loss (GN)0.02590.01010.04230.16020.1577 ScOT (50M) MSE loss 0.02330.00970.03960.14950.1581 Comp.loss (SA)0.02940.01090.05610.19520.1932 Comp.loss (GN)0.02890.01080.05970.19890.1992 DPOT-M (122M) MSE loss0.01950.00890.03450.12890.1396 Comp.loss (SA)0.01820.00820.03100.12090.1274 Comp.loss (GN) 0.01890.00860.03370.12490.1385 Poseidon-B (158M) MSE loss 0.07770.01760.49940.8943>10 Comp.loss (SA)0.07730.01670.52460.9230>10 Comp.loss (GN) 0.07890.01710.51110.9169>10 Table 33: Wavelet based spectral metrics for the 2D-SRBA (O) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining loss wfRMSE- high wfRMSE- mid wfRMSE- low MLW ConvNeXt (50M) MSE loss0.07090.16730.32750.8222 Comp.loss (SA) 0.06870.15920.31060.6784 Comp.loss (GN)0.06890.15840.30830.7178 CNO (50M) MSE loss 0.07920.19930.41410.7187 Comp.loss (SA)0.07210.17510.35110.5906 Comp.loss (GN)0.07540.18810.37940.6305 FFNO (50M) MSE loss0.06260.14610.29320.7409 Comp.loss (SA)0.05950.13840.27430.5820 Comp.loss (GN)0.06020.14020.27710.6177 ScOT (50M) MSE loss0.06450.14140.25650.7800 Comp.loss (SA)0.07170.16470.31310.7120 Comp.loss (GN)0.07140.16620.31380.7421 DPOT-M (122M) MSE loss0.05940.12430.21350.7759 Comp.loss (SA)0.05770.11930.20210.6703 Comp.loss (GN)0.06110.12650.21270.7185 Poseidon-B (158M) MSE loss0.10700.34090.93031.1640 Comp.loss (SA) 0.10360.33750.94401.0550 Comp.loss (GN)0.10570.34520.96141.0877 81 Table 34: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 2D-SRBA (O) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossIRMSESRMSE ConvNeXt (50M) MSE loss0.62970.2299 Comp.loss (SA)0.58210.2231 Comp.loss (GN) 0.57630.2230 CNO (50M) MSE loss0.74710.2975 Comp.loss (SA)0.62070.2582 Comp.loss (GN)0.65360.2795 FFNO (50M) MSE loss0.54460.2070 Comp.loss (SA)0.50780.1951 Comp.loss (GN) 0.51140.1975 ScOT (50M) MSE loss0.52050.1868 Comp.loss (SA)0.59560.2323 Comp.loss (GN) 0.59660.2345 DPOT-M (122M) MSE loss 0.15220.1571 Comp.loss (SA)0.14620.1554 Comp.loss (GN) 0.41030.1711 Poseidon-B (158M) MSE loss 0.84801.0545 Comp.loss (SA)0.83981.0778 Comp.loss (GN) 1.99871.0770 Table 35: Structure-aware metrics for the 2D-SRBA (O) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossPCCSSIM ConvNeXt (50M) MSE loss0.03000.1382 Comp.loss (SA) 0.02770.0806 Comp.loss (GN)0.02770.0891 CNO (50M) MSE loss0.04830.1911 Comp.loss (SA)0.03570.1047 Comp.loss (GN)0.04110.1262 FFNO (50M) MSE loss0.02690.0940 Comp.loss (SA)0.02500.0696 Comp.loss (GN)0.02520.0732 ScOT (50M) MSE loss 0.02060.0908 Comp.loss (SA)0.03050.0926 Comp.loss (GN)0.03160.1039 DPOT-M (122M) MSE loss0.01560.0854 Comp.loss (SA)0.01380.0518 Comp.loss (GN)0.01650.0640 Poseidon-B (158M) MSE loss0.45120.4469 Comp.loss (SA) 0.46840.4464 Comp.loss (GN)0.46590.4424 82 Table 36: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 2D-SRBA (O) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossiqRMSE-iqRMSE-iqRMSE-iqRMSE-iqRMSE- mass momentum _x momentum _y kinetic- energy enstrophy ConvNeXt (50M) MSE loss0.02770.0276>100.02660.0446 Comp.loss (SA)0.05410.0565>100.05210.0537 Comp.loss (GN)0.04200.0419>100.03760.0666 CNO (50M) MSE loss0.05660.0747>100.07730.0556 Comp.loss (SA) 0.03960.0546>100.05860.0551 Comp.loss (GN)0.04030.0517>100.05020.0657 FFNO (50M) MSE loss0.02940.0326>100.03450.0441 Comp.loss (SA)0.02470.0249>100.02320.0319 Comp.loss (GN)0.02550.0251>100.02300.0351 ScOT (50M) MSE loss 0.03690.0437>100.04680.0389 Comp.loss (SA)0.02340.0253>100.02500.0651 Comp.loss (GN) 0.02190.0234>100.02250.0874 DPOT-M (122M) MSE loss0.06240.0814>100.09510.0374 Comp.loss (SA)0.04970.0546>100.05280.0312 Comp.loss (GN)0.05060.0550>100.05430.0280 Poseidon-B (158M) MSE loss0.40953.3014>108.04321.6102 Comp.loss (SA)0.44253.1013>107.46231.4805 Comp.loss (GN) 0.37242.9242>107.33631.5545 Table 37: Pointwise evaluation of baselines on the 2D-SDBA (SSOO) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields and all rollout steps. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossH1H2MAERMSEVRMSE ConvNeXt (50M) MSE loss0.04110.01160.18910.47640.5366 Comp.loss (SA) 0.03800.01060.18540.46960.5242 Comp.loss (GN)0.03520.01000.17710.45520.5168 CNO (50M) MSE loss 0.04680.01340.21990.52070.6074 Comp.loss (SA)0.04030.01100.21110.50570.5755 Comp.loss (GN)0.03760.01070.19240.47360.5357 FFNO (50M) MSE loss0.03850.01060.18470.46150.5230 Comp.loss (SA)0.03750.01020.18220.46330.5204 Comp.loss (GN)0.03760.01030.17850.46240.5184 ScOT (50M) MSE loss0.04170.01380.17160.45130.5169 Comp.loss (SA)0.03490.01000.16410.43330.4969 Comp.loss (GN)0.03410.01000.16140.42470.4925 DPOT-M (122M) MSE loss0.04180.01270.14860.39620.4546 Comp.loss (SA)0.03170.00950.13230.36250.4061 Comp.loss (GN)0.03210.00950.13330.36560.4088 Poseidon-B (158M) MSE loss0.04530.01180.48470.81835.3877 Comp.loss (SA) 0.04350.01100.45430.78415.2695 Comp.loss (GN)0.04430.01050.49170.83275.8582 83 Table 38: Wavelet based spectral metrics for the 2D-SDBA (SSOO) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining losswfRMSE-wfRMSE-wfRMSE-MLW highmidlow ConvNeXt (50M) MSE loss0.08270.24030.64610.3711 Comp.loss (SA)0.07830.22770.61730.3456 Comp.loss (GN)0.07540.21620.58390.3465 CNO (50M) MSE loss0.08550.24880.68290.3773 Comp.loss (SA) 0.07920.23180.63590.3599 Comp.loss (GN)0.07680.22150.59720.3624 FFNO (50M) MSE loss0.08080.23500.62900.3657 Comp.loss (SA)0.07760.22710.61310.3385 Comp.loss (GN)0.07720.22630.61230.3291 ScOT (50M) MSE loss0.08750.23640.61680.3732 Comp.loss (SA) 0.07720.21800.58030.3400 Comp.loss (GN) 0.07620.21430.56380.3492 DPOT-M (122M) MSE loss0.08540.24080.60420.3599 Comp.loss (SA)0.07300.20520.52640.2963 Comp.loss (GN) 0.07320.20670.53350.3032 Poseidon-B (158M) MSE loss 0.08540.25660.72950.6353 Comp.loss (SA)0.08190.25010.71540.5974 Comp.loss (GN) 0.08070.25320.72950.6110 Table 39: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 2D-SDBA (SSOO) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossIRMSESRMSE ConvNeXt (50M) MSE loss 0.95500.4865 Comp.loss (SA)0.86390.4703 Comp.loss (GN) 0.79570.4505 CNO (50M) MSE loss 1.08930.5368 Comp.loss (SA)0.93060.5056 Comp.loss (GN)0.83880.4680 FFNO (50M) MSE loss 0.91700.4739 Comp.loss (SA)0.87190.4669 Comp.loss (GN)0.84560.4630 ScOT (50M) MSE loss0.91080.4606 Comp.loss (SA)0.78020.4315 Comp.loss (GN)0.77010.4264 DPOT-M (122M) MSE loss0.98220.4210 Comp.loss (SA)0.71480.3677 Comp.loss (GN)0.71470.3712 Poseidon-B (158M) MSE loss 1.33320.8327 Comp.loss (SA) 1.24490.7923 Comp.loss (GN)1.30530.8416 84 Table 40: Structure-aware metrics for the 2D-SDBA (SSOO) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossPCCSSIM ConvNeXt (50M) MSE loss0.17130.3158 Comp.loss (SA)0.16560.2949 Comp.loss (GN)0.16040.2846 CNO (50M) MSE loss0.20340.3572 Comp.loss (SA)0.18640.3205 Comp.loss (GN)0.16920.3043 FFNO (50M) MSE loss0.16900.3228 Comp.loss (SA) 0.16580.2957 Comp.loss (GN)0.16460.2874 ScOT (50M) MSE loss0.16450.3045 Comp.loss (SA)0.15140.2718 Comp.loss (GN)0.15150.2710 DPOT-M (122M) MSE loss 0.12630.2730 Comp.loss (SA)0.10400.2236 Comp.loss (GN)0.10470.2293 Poseidon-B (158M) MSE loss 0.33160.5316 Comp.loss (SA) 0.30470.4944 Comp.loss (GN)0.32180.5086 Table 41: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 2D-SDBA (SSOO) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossiqRMSE-iqRMSE-iqRMSE-iqRMSE-iqRMSE-iqRMSE- massmomentummomentumcom_xcom_yenstrophy _x_y ConvNeXt (50M) MSE loss 0.41690.21337.25690.03190.01800.3690 Comp.loss (SA)0.36530.22704.17690.01900.00950.3804 Comp.loss (GN) 0.28260.16573.79430.01570.01020.3865 CNO (50M) MSE loss>10>10>10>10>100.4177 Comp.loss (SA) 0.56450.6296>100.04880.02500.4205 Comp.loss (GN)>10>10>100.03900.04480.4259 FFNO (50M) MSE loss 0.19770.12193.40740.02300.01830.3737 Comp.loss (SA) 0.22080.18982.95570.02500.01630.3974 Comp.loss (GN)0.19210.19965.08660.01910.01210.3336 ScOT (50M) MSE loss >10>10>100.07750.08840.3506 Comp.loss (SA) 0.15830.10383.40000.01450.01090.3698 Comp.loss (GN)>10>10>100.03440.04070.4065 DPOT-M (122M) MSE loss >10>10>10>10>100.2815 Comp.loss (SA) >10>10>100.02140.02120.2072 Comp.loss (GN)>10>10>100.01600.01490.2134 Poseidon-B (158M) MSE loss0.07080.8916>100.03500.00983.3429 Comp.loss (SA)0.14650.8705>100.03150.00803.1110 Comp.loss (GN)0.04290.9015>100.03130.00503.5491 Table 42: Pointwise evaluation of baselines on the 3D-SABW (SSOOSS) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields and all rollout steps. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossH1H2MAERMSEVRMSE ConvNeXt MSE loss 0.05830.24280.14260.36971.9591 Comp.loss (SA) 0.06180.24340.17790.41901.6982 Comp.loss (GN)0.05340.22710.15760.37231.2033 FFNO MSE loss0.03330.18310.13180.30560.7823 Comp.loss (SA)0.03290.17390.15900.34930.6104 Comp.loss (GN)0.03030.16310.15270.33730.5767 85 Table 43: Wavelet based spectral metrics for the 3D-SABW (SSOOSS) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining losswfRMSE-highwfRMSE-midwfRMSE-lowMLW ConvNeXt MSE loss0.07290.33851.19392.1995 Comp.loss (SA)0.07900.34881.24372.0639 Comp.loss (GN)0.07340.31811.11021.9722 FFNO MSE loss0.06000.22200.76331.9471 Comp.loss (SA) 0.05700.22360.79821.7184 Comp.loss (GN)0.05500.21280.76171.7481 Table 44: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 3D-SABW (SSOOSS) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossIRMSESRMSE ConvNeXt MSE loss3.77610.4985 Comp.loss (SA)3.35040.5203 Comp.loss (GN)2.95740.4668 FFNO MSE loss1.79550.3214 Comp.loss (SA)1.60190.3427 Comp.loss (GN) 1.54210.3260 Table 45: Structure-aware metrics for the 3D-SABW (SSOOSS) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossPCCSSIM ConvNeXt MSE loss0.11600.4745 Comp.loss (SA)0.15890.4342 Comp.loss (GN) 0.11890.3721 FFNO MSE loss0.09230.3374 Comp.loss (SA)0.11070.3130 Comp.loss (GN) 0.10530.3165 Table 46: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 3D-SABW (SSOOSS) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossiqRMSE-iqRMSE-iqRMSE-iqRMSE-iqRMSE- mass momentum _x momentum _y kinetic- energy enstrophy ConvNeXt MSE loss0.28392.7939>10>100.1337 Comp.loss (SA)0.01700.0490>100.06820.2560 Comp.loss (GN) >10>10>10>100.1894 FFNO MSE loss0.00510.0218>100.03760.0508 Comp.loss (SA)0.00530.0187>100.02970.0673 Comp.loss (GN)0.00430.0118>100.01920.0546 86 Table 47: Pointwise evaluation of baselines on the 3D-SDBA (SSOOSS) test dataset, reporting errors in pointwise field values (MAE, RMSE, VRMSE) and the field derivatives (H1, H2) averaged across all fields and all rollout steps. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossH1H2MAERMSEVRMSE ConvNeXt MSE loss5.9460>100.29001.78214.3839 Comp.loss (SA)>10>100.24954.60786.0536 Comp.loss (GN)>10>10>10>10>10 FFNO MSE loss 0.11220.54870.15320.46051.0754 Comp.loss (SA)0.09590.43160.17610.46831.0765 Comp.loss (GN)0.09870.44910.18760.47441.4433 Table 48: Wavelet based spectral metrics for the 3D-SDBA (SSOOSS) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining losswfRMSE-highwfRMSE-midwfRMSE-lowMLW ConvNeXt MSE loss0.24531.20224.40121.3508 Comp.loss (SA)1.13133.21718.33661.0554 Comp.loss (GN) >10>10>101.1599 FFNO MSE loss0.10470.42031.39611.1078 Comp.loss (SA)0.08880.37411.31791.0079 Comp.loss (GN)0.09030.38581.32991.0078 Table 49: Feature focused metrics for capturing interface (IRMSE) and shocks (SRMSE) on the 3D-SDBA (SSOOSS) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossIRMSESRMSE ConvNeXt MSE loss >102.0698 Comp.loss (SA) >103.8508 Comp.loss (GN)>10>10 FFNO MSE loss2.05070.5131 Comp.loss (SA) 1.34020.4788 Comp.loss (GN)1.56610.5041 Table 50: Structure-aware metrics for the 3D-SDBA (SSOOSS) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossPCCSSIM ConvNeXt MSE loss0.34390.5562 Comp.loss (SA)0.24860.3653 Comp.loss (GN)0.31240.4845 FFNO MSE loss 0.16610.3790 Comp.loss (SA)0.16120.3620 Comp.loss (GN)0.16910.3794 87 Table 51: Physics-informed metrics reporting the nRMSE of integral quantities of interest on the 3D-SDBA (SSOOSS) test dataset. Lower values are better. Darker shade indicates lower value for each baseline within a group of training loss strategies. Bold values denote the lowest among all baselines for a given metric and the underlined values denote the second-lowest values. ModelTraining lossiqRMSE-iqRMSE-iqRMSE-iqRMSE-iqRMSE-iqRMSE- massmomentummomentumcom_xcom_yenstrophy _x_y ConvNeXt MSE loss>10>10>10>10>100.9275 Comp.loss (SA)>10>10>10>10>100.8792 Comp.loss (GN)>10>10>10>10>101.1215 FFNO MSE loss>10>10>100.09850.08280.2710 Comp.loss (SA) >10>10>100.05220.03580.3176 Comp.loss (GN)2.41732.2959>100.04870.03000.3190 88 C.4 Rollout visualization Density Prediction TargetAbs. Error Input Pressure Prediction TargetAbs. Error Input 2D-SABW-O ConvNeXt-50M Figure 69: Rollout visualization of Density and Pressure fields for 2D Shock-induced Air Bubble Collapse in Water with Open Boundaries (2D-SABW-O) using ConvNeXt-50M. 89 VelocityX Prediction TargetAbs. Error Input VelocityY Prediction TargetAbs. Error Input 2D-SABW-O ConvNeXt-50M Figure 70: Rollout visualization of Velocity-X and Velocity-Y fields for 2D Shock-induced Air Bubble Collapse in Water with Open Boundaries (2D-SABW-O) using ConvNeXt-50M. 90 Density Prediction TargetAbs. Error Input Pressure Prediction TargetAbs. Error Input 2D-SRBA-O CNO-50M Figure 71: Rollout visualization of Density and Pressure fields for 2D Shock-induced R22 Bubble Collapse in Air with Open Boundaries (2D-SRBA-O) using CNO-50M. 91 VelocityX Prediction TargetAbs. Error Input VelocityY Prediction TargetAbs. Error Input 2D-SRBA-O CNO-50M Figure 72: Rollout visualization of Velocity-X and Velocity-Y fields for 2D Shock-induced R22 Bubble Collapse in Air with Open Boundaries (2D-SRBA-O) using CNO-50M. 92 Density Prediction TargetAbs. Error Input Pressure Prediction TargetAbs. Error Input 2D-SABW-O FFNO-50M Figure 73: Rollout visualization of Density and Pressure fields for 2D Shock-induced Air Bubble Collapse in Water with Open Boundaries (2D-SABW-O) using FFNO-50M. 93 VelocityX Prediction TargetAbs. Error Input VelocityY Prediction TargetAbs. Error Input 2D-SABW-O FFNO-50M Figure 74: Rollout visualization of Velocity-X and Velocity-Y fields for 2D Shock-induced Air Bubble Collapse in Water with Open Boundaries (2D-SABW-O) using FFNO-50M. 94 Density Prediction TargetAbs. Error Input Pressure Prediction TargetAbs. Error Input 2D-SDBA-SSOO-SIE ScOT-50M Figure 75: Rollout visualization of Density and Pressure fields for 2D Shock-induced Droplet Breakup in SIE regime with Symmetric Boundaries (2D-SDBA-SSOO) using ScOT-50M. 95 VelocityX Prediction TargetAbs. Error Input VelocityY Prediction TargetAbs. Error Input 2D-SDBA-SSOO-SIE ScOT-50M Figure 76: Rollout visualization of Velocity-X and Velocity-Y fields for 2D Shock-induced Droplet Breakup in SIE regime with Symmetric Boundaries (2D-SDBA-SSOO) using ScOT-50M. 96 Density Prediction TargetAbs. Error Input Pressure Prediction TargetAbs. Error Input 2D-SDBA-SSOO-RTP DPOT-M-122M ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... Figure 77: Rollout visualization of Density and Pressure fields for 2D Shock-induced Droplet Breakup in RTP regime with Symmetric Boundaries (2D-SDBA-SSOO) using DPOT-M-122M. 97 VelocityX Prediction TargetAbs. Error Input VelocityY Prediction TargetAbs. Error Input 2D-SDBA-SSOO-RTP DPOT-M-122M ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... Figure 78: Rollout visualization of Velocity-X and Velocity-Y fields for 2D Shock-induced Droplet Breakup in RTP regime with Symmetric Boundaries (2D-SDBA-SSOO) using DPOT-M-122M. 98 Density Prediction TargetAbs. Error Input Pressure Prediction TargetAbs. Error Input 2D-SABW-SSOO Poseidon-B-158M Figure 79: Rollout visualization of Density and Pressure fields for 2D Shock-induced Air Bubble Collapse in Water with Symmetric Boundaries (2D-SABW-SSOO) using Poseidon-B-158M. 99 VelocityX Prediction TargetAbs. Error Input VelocityY Prediction TargetAbs. Error Input 2D-SABW-SSOO Poseidon-B-158M Figure 80: Rollout visualization of Velocity-X and Velocity-Y fields for 2D Shock-induced Air Bubble Collapse in Water with Symmetric Boundaries (2D-SABW-SSOO) using Poseidon-B-158M. 100 NeurIPS Paper Checklist 1. Claims Question: Do the main claims made in the abstract and introduction accurately reflect the paper’s contributions and scope? Answer: [Yes] Justification: The claims made in the abstract are covered in sections 3 and 5 of the main text. Further information are detailed in the appendix section C. 2. Limitations Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: We have provided the limitations clearly in the section 5.2 of the main text. 3. Theory assumptions and proofs Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [N/A] Justification: The paper does not include theoretical results, but rather focuses on the empirical benchmarking of surrogate models for shock-induced multiphase flows. 4. Experimental result reproducibility Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: For generating the datasets, we provide all the details related to the physics and the solver in the metadata.json file. For reproducing the benchmarks we provide the model hyperparameters in section B.1, training hyperparameters in section B.2 and hyperparameters related to metrics and loss-weighting strategies in section B.4. Further the hardware specifications on which the training was carried out is specified in section B.6. 5. Open access to data and code Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: The link to the open-sourced (anonymized) code and datasets is provided in the footnote of the first page of the main text. The README.md of the benchmarking repository contains all the necessary instructions to set up a virtual environment, train the baselines and generate metrics from the trained models during inference. 6. Experimental setting/details Question: Does the paper specify all the training and test details (e.g., data splits, hyperpa- rameters, how they were chosen, type of optimizer) necessary to understand the results? Answer: [Yes] Justification: The data-splits are provided in Figure 1, the hyperparameters for the baselines in section B.1, training hyperparameters in section B.2 and hyperparameters related to metrics and loss-weighting strategies in section B.4 7. Experiment statistical significance Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] Justification: We provide error bars in our metric evolution plots for the 2D and 3D baselines in section C.2. 101 8. Experiments compute resources Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer:[Yes], Justification: The wall clock time required to generate the datasets are specified in Table 2. The hardware specifications for generating the datasets and the models is specified in section B.6. The peak memory requirements and the wall clock time for a single training epoch is presented in C.1. 9. Code of ethics Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.c/public/EthicsGuidelines? Answer: [Yes] Justification: We reviewed the NeurIPS Code of Ethics, and we confirm that our research conforms to it in every respect. 10. Broader impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [N/A] Justification: The datasets and trained models are restricted to simplified benchmark con- figurations and are not directly deployable in real-world decision-making. As such, no immediate societal risks or negative impacts are anticipated, and the broader societal impact of this work is expected to remain primarily within the research community. 11. Safeguards Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pre-trained language models, image generators, or scraped datasets)? Answer: [N/A] Justification: The datasets and the trained surrogates provided in this paper are intended for scientific benchmarking and methodological evaluation of surrogate models rather than deployment in safety-critical or real-world operational systems 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: The solver used to generate the datasets is open sourced and is cited in the paper and a link to the repository-https://gitlab.lrz.de/nanoshock/ALPACA. The license of the solver is GNU General Public License v3.0. 13. New assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [Yes] Justification: The datasets and benchmarking pipeline were developed by the authors. The datasets are publicly hosted on Hugging Face, while the benchmarking repository is maintained on GitHub. All required licensing information is clearly specified in the README.md files. The baseline models used for benchmarking are adapted from existing open-source implementations, and appropriate credits and citations are provided within the benchmarking repository. 14. Crowdsourcing and research with human subjects 102 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)? Answer: [N/A] Justification: The paper does not involve crowdsourcing nor research with human subjects. 15.Institutional review board (IRB) approvals or equivalent for research with human subjects Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? Answer: [N/A] Justification: The paper does not involve crowdsourcing nor research with human subjects. 16. Declaration of LLM usage Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigor, or originality of the research, declaration is not required. Answer: [N/A] Justification: LLMs were used only for writing, editing, and formatting purposes and did not impact the core methodology, scientific rigor, or originality of the research. 103