Spectral-Embedded Operator Learning for Three-Phase Interfacial Flow: A Ternary Cahn-Hilliard-Navier-Stokes Benchmark

Muhammad Abid; Arth Sojitra; Omer San

Source-linked AI summary

Spectral-Embedded Operator Learning for Three-Phase Interfacial Flow: A Ternary Cahn-Hilliard-Navier-Stokes Benchmark

Muhammad Abid, Arth Sojitra, Omer San

arXiv:2608.29069v1Aug 29, 2026physics.flu-dyncs.LG

TL;DR

Existing operator-learning evaluations leave unclear whether trunk coordinate choices transfer to constrained, multiphase interface-dominated flows. This paper introduces a structure-preserving ternary benchmark and compares matched DeepONet trunks, finding that Chebyshev embeddings outperform raw and Fourier coordinates on the tested non-periodic flow.

Problem
Existing spectral-trunk evidence largely concerns single-field targets on unconstrained domains, leaving transfer to constrained three-phase interface-dominated fields unresolved.
Method
The paper learns a nine-parameter-to-five-channel operator from 1,024 simulations using a structure-preserving ternary solver and parameter-matched DeepONet variants.
Results
24.0% lower test relative L2 error than DeepONet and 16.8% lower than FEDONet, with improvements across all five output channels.
Takeaways & Limitations
Chebyshev trunk coordinates are particularly effective for the tested wall-confined, strongly non-periodic flow, with gains concentrated near interfaces and after breakthrough.
Takeaways & Limitations
A rank-limited representation remains an accuracy bottleneck for translating sharp interfaces, and coordinate dictionaries do not change the reachable basis.

Abstract

from arXiv · show

Operator-learning surrogates have been benchmarked largely on single-field, single-interface problems, leaving unclear whether architectural choices validated in those settings transfer to constrained, multiphase flows. We introduce a three-phase interfacial-flow benchmark to examine whether the trunk coordinate representation matters for a multi-channel, interface-dominated target. The configuration consists of an air bubble rising through water, piercing a water-oil interface, and entraining a water plume into the oil within a bounded, wall-confined domain. Reference data are generated using a structure-preserving ternary Cahn-Hilliard-Navier-Stokes solver that algebraically preserves the simplex constraint. From 1,024 Sobol-sampled simulations spanning a nine-dimensional parameter space, we learn the mapping from physical parameters to five-channel space-time fields. We compare three parameter-matched DeepONet variants differing only in trunk representation: raw coordinates (DeepONet), random Fourier features (FEDONet), and a fixed tensor-product Chebyshev dictionary (SEDONet). SEDONet reduces the test relative L2 error by 16.8% compared with FEDONet and by 24.0% compared with DeepONet, while improving all five output channels. Spatial and temporal error analyses localize the principal gains near the diffuse interfaces and after bubble breakthrough. The results indicate that the Chebyshev representation is particularly effective for the strongly non-periodic wall-normal and temporal structure of this three-phase flow.

1 Introduction

Three-phase interfacial flow combines multiple coupled interfaces, topological transitions, algebraic phase constraints, and expensive simulation, creating a demanding test for operator-learning architectures. This paper introduces a controlled benchmark and compares trunk coordinate representations on that target.

Three immiscible fluids create three pairwise interfaces, jointly determined triple-junction geometry, and regime boundaries from partial or total spreading.
The benchmark targets an air bubble that rises through water, crosses a water-oil interface, and entrains a thinning, pinching water plume.
Diffuse-interface simulation resolves topology changes automatically but requires fine interface resolution and restrictive time steps, making parametric studies expensive.
Existing spectral-trunk evidence largely concerns single-field, unconstrained problems, leaving transfer to constrained five-channel three-phase targets unresolved.
The benchmark contains a five-channel space-time field over a nine-dimensional parameter space generated from 1,024 Sobol-sampled simulations.
SEDONet, FEDONet, and DeepONet are compared with matched training components, while SEDONet uses a degree-balanced Chebyshev dictionary for three coordinates.

2 Related Work

Prior work establishes operator learning and spectral trunk embeddings but has mostly evaluated simpler, unconstrained, or periodic settings. The paper positions its benchmark as a test of these choices under three-phase interfaces, non-periodic boundaries, and structural constraints.

Sharp-interface methods achieve strong two-phase geometric fidelity but require dedicated junction models and topology-changing surgery in three-phase settings.
Diffuse-interface methods make merging and pinching automatic, but their interface width must be resolved by the computational grid.
Prior operator-surrogate successes largely involve single-phase or effectively two-phase, periodic or quasi-periodic, unconstrained targets.
Neural spectral bias motivates coordinate embeddings, while random Fourier features reshape the learned spectrum and have shown gains on oscillatory PDEs.
Chebyshev bases are suited to bounded non-periodic intervals because they converge geometrically for analytic functions and cluster nodes near endpoints.
The benchmark preserves the simplex constraint in reference data so surrogate violations can be evaluated as model errors rather than data inconsistencies.
A translating sharp feature remains difficult for every fixed-basis architecture, limiting what coordinate dictionaries can ultimately represent.

3 Governing Model: Ternary Cahn–Hilliard–Navier–Stokes

The reference model is a bounded-domain ternary Cahn–Hilliard–Navier–Stokes system for water, oil, and air, with diffuse interfaces, variable mixture properties, and an exact simplex constraint. Mobility choices and homogeneous Neumann conditions make that constraint and transform-based discretization structural features of the solver.

The model produces a bubble rising through water, crossing a water-oil interface, and entraining a plume in a bounded wall-confined domain.
Water, oil, and air volume fractions lie on the Gibbs simplex pointwise, making the three phase fields algebraically dependent.
Positive spreading coefficients define partial spreading with triple-junction angles jointly fixed by the three surface tensions, whereas negative coefficients produce total spreading.
The free-energy formulation uses a gradient-free chemical-potential component and a mobility choice that makes the weighted chemical-potential contribution cancel identically.
Variable-density and variable-viscosity mixture properties couple phase composition to Navier–Stokes dynamics, with buoyancy driving bubble rise.
The potential-form interfacial force avoids explicit curvature and triple-junction modeling, which helps retain tractability through topology changes.
Homogeneous Neumann scalar conditions and no-slip velocity conditions support exact cosine diagonalization of the discrete Laplacian in the strictly non-periodic domain.
Choosing Mi = M0/Σi propagates the simplex constraint exactly through the discretized system rather than enforcing it with penalties or projection.

4 Operator-Learning Formulation

The formulation maps nine physical parameters to constrained five-channel space-time fields using a shared branch–trunk basis, while varying only the fixed trunk coordinate embedding. Its comparison is bounded by rank-p representation limits and tests raw, Fourier, and Chebyshev coordinate dictionaries under controlled conditions.

4.1 The parametric solution operator: The learned operator maps a nine-dimensional parameter vector to an entire five-channel space-time field, rather than predicting one discretized solution.The branch consumes the nine parameters directly, so sensor-resolution error is absent; measured error is representation plus optimization error.
4.1 The parametric solution operator: The target operator is required to be single-valued and continuous across the sampled parameter box, excluding ensemble members with qualitatively different topological outcomes.Diagnostics rule out cases such as failed bubble crossing or alternative plume behavior within the ensemble.
4.2 Branch-Trunk factorization and the simplex closure: The rank-p factorization imposes a Kolmogorov p-width error floor that no trunk embedding or additional training can overcome.Moving interfaces are poorly represented by fixed linear modes, motivating nonlinear basis parameterizations as future work.
4.2 Branch-Trunk factorization and the simplex closure: The shared trunk supplies spatiotemporal basis functions for all five channels, while the branch produces channel-specific coefficients and the simplex closure reweights errors affecting the air-interface location.Sharing encodes common interface geometry without requiring each channel to rediscover it independently.
4.3 Trunk embeddings: the only difference between the three models: All three models retain identical networks, losses, optimization, and data, differing only in the fixed embedding of affinely rescaled query coordinates.Rescaling makes the spatial and temporal axes comparable despite their different physical extents.
4.3.2 FEDONet: Random Fourier features: Random Fourier features impose a stationary, isotropic, periodic prior that is mismatched to bounded non-periodic structure and endpoint behavior.The representation treats walls like interior points and can require slowly converging sinusoidal expansions for linear trends.
4.3.3 SEDONet: Chebyshev dictionary: The Chebyshev trunk uses an orthogonal polynomial dictionary on the rescaled interval and orders tensor-product modes by total degree to balance coordinate resolution.This avoids the lexicographic crop’s near-blindness to higher-order structure in the first coordinate.
4.4 Why the embedding changes the optimization: Chebyshev features represent linear trends exactly in two modes, unlike random Fourier features, which incur endpoint Gibbs penalties, or raw nonlinear trunks.The distinction targets monotone wall-normal and temporal behavior such as bubble rise, plume advance, and entrained-volume growth.

5 Results

SEDONet improves test accuracy across all five channels, with gains that are broad across the design space and concentrated at diffuse interfaces after bubble breakthrough. The results support Chebyshev features for monotone, non-periodic wall-normal and temporal structure while showing that Fourier features retain advantages for some periodic, oscillatory structures.

Test accuracy: 24.0% lower mean relative L2 error than the raw-trunk baseline and 16.8% lower than FEDONet, uniformly across all five channels.The shared trunk exposes whether a representation helps some channels at the expense of others; SEDONet does not show that trade-off.
Per-channel accuracy: SEDONet improves the two large-region concentration fields cw and co by about a fifth while also exceeding Fourier gains on oscillatory channels.FEDONet’s gains concentrate mainly on ca and the velocities, whereas SEDONet improves every channel.
Mechanism: The polynomial dictionary represents monotone, strongly non-periodic trends efficiently, leaving more capacity for interface structure than a periodic sinusoidal dictionary.A linear ramp requires only two Chebyshev modes but a slowly converging sinusoidal series under the periodic representation.
Across the design: SEDONet remains below the other models across the three most influential parameters, with roughly constant logarithmic separation rather than gains confined to a design subregion.The improvement is therefore broad across the sampled parameter space, although individual-member errors still vary substantially.
Where the error lives: The principal spatial gains occur at moving interface peaks, where the baseline reaches 0.33 at t = 4 versus 0.17 for SEDONet, while bulk errors remain near zero.The models separate at the bubble rim, deformed interface, and plume flanks rather than throughout the liquid bulk.
When the error emerges: The models are similar before breakthrough, but their temporal ordering reverses at breakthrough and SEDONet’s advantage widens thereafter.Before breakthrough the target is smooth and effectively low-rank; afterward it contains a triple junction, a thin high-curvature spike, and pronounced monotone drift.
Hard test member: On the hard test member, SEDONet is lowest on every channel by 37.0% versus the baseline and 19.3% versus FEDONet at t = 2.At t = 4, the corresponding margins are 29.3% and 10.6%; at t = 3, FEDONet briefly wins on ca while SEDONet retains the other listed channels.
Qualitative fields: All models reproduce the correct topology, and their differences appear mainly as the thickness and brightness of a thin error shell around interfaces.The remaining error is characterized as interface-localization error rather than a bulk field-level failure.

6 Summary and Conclusions

The benchmark compares parameter-matched DeepONet variants on a constrained, five-channel ternary-flow target. SEDONet achieves the strongest overall results, while the evidence supports a narrower interpretation tied to non-periodic coordinate structure and leaves important validation limits.

The benchmark adds a five-channel space-time target, an algebraic state constraint, and a structurally constraint-preserving reference solver.
24.0% lower test relative L2 error than DeepONet and 16.8% lower than FEDONet were obtained by SEDONet.The comparison was parameter-exact: the models shared depths, widths, optimizers, data, and parameter counts, differing only in trunk embedding.
SEDONet improved every output channel and the bubble-centroid diagnostic while preserving the simplex constraint to the single-precision rounding floor.
The advantage is localized to diffuse interfaces and post-breakthrough behavior rather than being attributed solely to boundary resolution.Errors are smallest at the walls, while the proposed interpretation concerns monotone, non-periodic wall-normal and temporal trends.
The reported gaps are point estimates from one training run, and transfer beyond finite-Cahn-number diffuse solutions or one topological class remains limited.

7 Future Work

Future work prioritizes stronger evaluation and broader physical coverage, alongside architectures that adapt to moving interface geometry. The proposed extensions target uncertainty in current effect sizes and the limits of fixed global dictionaries.

Repeated runs at matched compute budgets, multiple seeds, and varied Fourier bandwidths are needed to make the reported effect sizes more defensible.
The design space should cross regime boundaries and extend beyond the training box to test generalization across topological and parameter changes.
A dictionary combining global trend modes with localized modes tied to instantaneous interface position could address the error concentrated at diffuse interfaces.

A Proof of Proposition 3.1

The proof establishes exact propagation of the simplex constraint by cancellation in the continuous and compatible discrete formulations. Its conditioning appendix also reports exact weighted-grid behavior for the Chebyshev dictionary.

Summing the phase equations with the prescribed weights cancels the relevant terms and yields exact propagation of Σ_i c_i ≡ 1.
The discrete cancellation remains valid when the same discrete Laplacian is used in the coupled equations.
Chebyshev tensor-product modes are orthogonal under the stated weighted inner product, with normalization constants γ_0 = π and γ_n≥1 = π/2.
The discrete Gauss–Chebyshev counterpart is the property seen by the embedding on the tensor-product grid.
With K = 10, the weighted-grid Gram matrix has off-diagonal entries at most 1.8 × 10^-16 and condition number 2 per axis, giving 8 in three dimensions.

B.2 Conditioning on the uniform query grid

On the uniformly sampled query grid, Chebyshev conditioning is weaker than under its orthogonality weight but remains deterministic and fixed. The comparison with Legendre does not establish a polynomial-family advantage.

Uniform query sampling produces per-axis condition numbers of 15.1, 13.7, and 5.2 for the evaluated Chebyshev dictionary.
Grid refinement drives the per-axis condition figure toward 13.51 rather than the weighted-grid value 2.
The realized whitening is weaker by about thirtyfold, but the Gram matrix remains fixed, deterministic, seed-identical, and bandwidth-free.
Nothing in the mechanism selects Chebyshev over Legendre, and that comparison remains untested.
The Fourier comparison is described as comparative and requires the same conditioning treatment.

C.1 Induced kernel

The Fourier-feature trunk induces a stationary, isotropic prior, while the Chebyshev dictionary has boundary-enhanced structure and fixed conditioning. Fourier conditioning varies sharply with bandwidth, including collapse at small σB and improvement near σB = 4.

Fourier-feature prior: The Fourier kernel is stationary and isotropic, imposing the same prior at walls and interior and one resolvable scale across all three axes.The polynomial prior is maximal at the corners of [−1, 1]3, whereas the Fourier prior does not distinguish boundary from interior.
Conditioning: κ(G) ≈2.7 × 102 for the graded Chebyshev dictionary, with no parameter to select.
Fourier-feature prior: Distinct frequency rows are uncorrelated and yield E[G] = 1 2I, so the Fourier dictionary is whitened in expectation.The optimizer nevertheless sees the realized Gram matrix, whose conditioning depends sharply on σB.
Conditioning: As σB →0, Fourier features approach rank one and G becomes singular to single precision.
Conditioning: Near σB = 4, the Fourier dictionary is better conditioned than the polynomial dictionary by two orders of magnitude.

C.3 The cost of a monotone trend

Chebyshev modes represent monotone wall-normal and temporal trends far more efficiently than Fourier features at the same budget. Increasing Fourier bandwidth improves conditioning but worsens trend approximation, creating a structural trade-off within that family.

Trend representation: 6 × 10−17 relative residual: the first two Chebyshev modes fit f(ξ) = ξ exactly on the wall-normal grid.This follows from T1 = ξ.
Trend representation: 6.9 × 10−3, 2.8 × 10−2, 1.2 × 10−1 and 6.3 × 10−1 are the Fourier residuals at σB = 0.5, 1, 2 and 4, respectively.
Trend representation: 8 features at σB = 1 and 16 at σB = 2 are required for Fourier residuals below 10−3.
Trade-off: Fourier residuals grow with bandwidth while Gram-matrix conditioning improves, so no single σB serves both objectives.The trade-off is structural rather than merely a tuning issue; the trend mismatch is not claimed to explain an order-of-magnitude gap.

E The benchmark dataset

The benchmark uses a structure-preserving, fixed-cost numerical pipeline designed to generate 1,024 simulations efficiently while retaining discrete simplex structure. Boundary-aware transforms and species-independent operators remove iterative solves and preserve the algebraic closure.

Benchmark specification: Tables 9 and 10 together specify the benchmark discretization, fixed physical quantities, cost, and nine-variable design.
Discretization: The reference solver uses staggered MAC finite differences, second-order centered space, first-order time, explicit nonlinear terms, and an implicit linear Cahn–Hilliard part.
Efficient solver: Homogeneous Neumann conditions make the discrete Laplacian exactly DCT-II diagonalizable, turning the pressure and linear Cahn–Hilliard solves into one-shot operations.These solves require no iteration, tolerance, or run-to-run cost variability.
Simplex structure: MiΣi = M0 and MiSi = M0S0 make the implicit operator species-independent across all three phases.This common left-hand operator carries the continuous simplex-preservation result to the discrete level.
Simplex structure: The discrete simplex proposition assumes divergence-free velocity, linear face interpolation, and pointwise phase closure P_i = 1.Summing the species equations then yields the closure pointwise and unconditionally.
Ensemble execution: One jitted JAX kernel is vmapped over all 1,024 members, with one compilation completing the ensemble in 37 minutes.The pressure correction must reuse the lagged pressure from the Poisson right-hand side to preserve the stated discrete properties.

E.2 Initial condition and ensemble design

Each simulation starts from a resting circular air bubble below a flat interface, with the third phase obtained by affine closure. The ensemble contains 1,024 admissible Sobol designs spanning nine variables and uses a domain sized to contain the plume.

Initial condition: Each member starts at rest with a circular air bubble of radius d/2 centered below a flat interface at yint.
Initial condition: The third phase is constructed by affine closure, so the simplex theorem applies from the first time step.Because the initial surrogate map is a closed-form function of π, its nonzero t = 0 error is attributed to trunk representation.
Ensemble design: 1,024 = 2^10 scrambled Sobol points span nine parameters, with no rejected draws because Σi > 0 throughout the sampling box.The realized design therefore retains the sequence's low-discrepancy construction.
Ensemble design: Parameters 1–3 control capillarity and wetting, parameters 4–7 material contrasts, and parameters 8–9 initial bubble placement.

E.3 Diagnostics, well-posedness and known imperfections

The diagnostics establish a well-posed, single-topological-class operator-learning target while documenting coarse-grid and physical-model imperfections. Across the ensemble, cases follow a consistent rise-to-detachment sequence and exhibit monotone responses over the Eo sweep.

Diagnostics: Four scalar diagnostics are stored at every snapshot to verify that the ensemble is well posed and provide Section 5.4 accuracy measures.
Ensemble behavior: Every ensemble member follows the same qualitative sequence: rise, deform, breakthrough, entrain, and detach.
Well-posedness: A sweep across Eo ∈[5, 40] produces a monotone response, with every case crossing the interface and entraining a plume.
Well-posedness: The single topological class makes G single-valued, but the surrogates are never tested near a regime boundary.
Known imperfections: Coarsening preserves the simplex constraint exactly, whereas stored velocities remain divergence free only to interpolation accuracy.
Known imperfections: The ensemble exhibits at most 2.1 × 10−4 simplex drift and loses 15% of its initial bubble area on average by t = 4.