cs.LGOct 21, 2026

Local and Sparse Kernel Operator Learning for Incompressible Flows: An Operator Lagrange Formulation

Authors: Ramansh Sharma, Matthew Lowery, Houman Owhadi, Varun Shankar

Organizations: Kahlert School of Computing, University of Utah, UT, USA · Department of Computing and Mathematical Sciences, California Institute of Technology, CA, USA

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

We present a local and sparse kernel method for learning the solution operators of incompressible flows. It builds on our earlier property-preserving kernel method, whose global output reconstruction and dense operator regression scale poorly. Our method replaces the global output reconstruction by local divergence-free interpolants with divergence-free polynomials of any degree, solved in parallel on the GPU. These interpolants are divergence free, periodic on periodic domains, and, with a multiscale kernel, approach Kolmogorov statistics under refinement. A compactly supported kernel, positive semidefinite in any input dimension, makes the operator regression sparse. Its support radius, and hence the sparsity, is fitted to data, and optional affine polynomials in the input capture global trends. Operator Lagrange weights join the two. They cost one sparse solve per new input and act on the training fields before the output reconstruction, so their local coefficients are never formed. We show that, unlike direct evaluation, Lagrange evaluation does not amplify solve errors by the cancellation among spatial coefficients. On ten fluid benchmarks, validation selects sparse matrices nearly everywhere, with densities that mostly decrease as the training set grows, and affine augmentation lowers the error by more than an order of magnitude for forced isotropic turbulence. These sparse models are at least as accurate as our earlier dense method and more accurate than neural operators. Lagrange evaluation is as accurate as direct evaluation with local output reconstructions, and up to nearly three orders of magnitude more accurate with global ones. We also present GPU timings and storage.

Figures & tables

Figure 1: Overview of our method for a new input a∗\mathbf{a}_{*} . The compactly supported kernel λ\lambda of Section 3 gives a sparse matrix KUUθK_{UU}^{\theta} on the sampled training inputs (left and center). One sparse solve with the kernel values kU(a∗)k_{U}(\mathbf{a}_{*}) gives the operator Lagrange weights ℓ(a∗)\ell(\mathbf{a}_{*}) of Section 4 , shown here without their affine correction. The same weights give the coefficients b∗j=BjTℓ(a∗)\mathbf{b}_{*}^{j}=B_{j}^{\mathsf{T}}\ell(\mathbf{a}_{*}) of the local divergence-free interpolant of Section 2 on every stencil SjS_{j} of the output (right), where the rows of BjB_{j} are the coefficients of the training fields on SjS_{j} . In practice, we apply the weights to the sampled training fields before the output reconstruction, which gives the same prediction (Section 4.2 ).
StepWorkStorage
Kernel factorizationFcholF_{\rm chol}NK\mathcal{N}_{K} , Nfac\mathcal{N}_{\rm fac}
Affine correction ( naffn_{\rm aff} columns)solves with naffn_{\rm aff} right-hand sides, O(Nnaff2)O(Nn_{\rm aff}^{2}) QRO(Nnaff)O(Nn_{\rm aff})
Lagrange weights, bqb_{q} queriesO(bqNfac)+O(bqNnaff)O(b_{q}\mathcal{N}_{\rm fac})+O(b_{q}Nn_{\rm aff})bqNb_{q}N (dense)
Combining training fieldsO(bqNnv)O(b_{q}Nn_{v})NnvNn_{v}
Local systems (setup)O(rj3)O(r_{j}^{3}) per stencil 
Local evaluation, m⋆m^{\star} pointsO(bqm⋆d2ns)O(b_{q}m^{\star}d^{2}n_{s}) per snapshotO(m⋆d2ns)O(m^{\star}d^{2}n_{s}) weights, O(m⋆ns)O(m^{\star}n_{s}) indices
Table 1: Work and storage. Here NK\mathcal{N}_{K} and Nfac\mathcal{N}_{\rm fac} are the numbers of nonzeros of the kernel matrix and of its Cholesky factor, FcholF_{\rm chol} the work ( 68 ) of the factorization, bqb_{q} the number of new inputs in a batch, naffn_{\rm aff} the number of affine columns, nvn_{v} the number of sampled output values per field, m⋆m^{\star} the number of query points, rjr_{j} the size of the local system on stencil jj , and mrhs=∑jrjm_{\rm rhs}=\sum_{j}r_{j} . The last two rows belong to direct evaluation through coefficients, which our predictions avoid. For spacetime outputs, the counts of these two rows are multiplied by NtN_{t} .
Figure 3: Convergence of the local output reconstruction on scattered nodes in the domains Ω2\Omega_{2} (left) and Ω3\Omega_{3} (right), with PHS and multiscale Matérn kernels and divergence-free polynomials. Relative error EspE_{\rm sp} at random evaluation points that are not nodes, against the sampled fill distance hQh_{Q} of ( 69 ). Dotted lines: hQ2h_{Q}^{2} , hQ4h_{Q}^{4} , hQ6h_{Q}^{6} .
Figure 4: Kolmogorov statistics of the local multiscale Matérn output reconstruction. Left: structure functions for degree five on three of the six node sets, and the exact reference. Right: worst relative error ES2(0.8)E_{S_{2}}(0.8) of ( 71 ) of the reconstructed structure function against the node spacing hh , with an h2/3h^{2/3} guide.
Data setNodes mmEvaluation pointsGlobalLocal
Forced isotropic turbulence700030001.88×10−21.88\times 10^{-2}1.65×10−2\mathbf{1.65\times 10^{-2}} (multiscale)
Turbulent species transport700043881.52×10−11.52\times 10^{-1}1.47×10−1\mathbf{1.47\times 10^{-1}} (multiscale)
Merging vortices50069471.88×10−11.88\times 10^{-1}1.35×10−1\mathbf{1.35\times 10^{-1}}
Lid-driven cavity100085666.15×10−16.15\times 10^{-1}8.71×10−2\mathbf{8.71\times 10^{-2}}
Cylinder flow (no shedding)100085208.67×10−2\mathbf{8.67\times 10^{-2}}9.19×10−29.19\times 10^{-2}
Cylinder flow (vortex shedding)100085209.65×10−2\mathbf{9.65\times 10^{-2}}1.12×10−11.12\times 10^{-1}
Table 2: Test error between the nodes. We reconstruct the predicted nodal velocities of the 100 test inputs (model with the smallest validation error at N=10000N=10000 ) and measure the mean relative velocity error at the evaluation points, the mesh points of each data set that are not among its mm nodes. Global: the divergence-free Matérn interpolant of Sharma et al. [53] on all nodes. Local: our output reconstruction with the cubic PHS, or with the multiscale Matérn kernel ( cσ=16c_{\sigma}=16 ) where marked (Appendix A.5 compares it with the PHS). Bold: smaller error in each row.
Figure 5: Setup time and memory of our local output reconstruction (degree four, cubic PHS) and of a global divergence-free Matérn output reconstruction (dense Gram matrix and Cholesky factor) against the number mm of output points, on scattered nodes in the unit square (2D) and cube (3D). Dotted: the global cost extrapolated as (dm)3(dm)^{3} in time and (dm)2(dm)^{2} in memory beyond the largest size that a single dense factorization fits on our GPU. All times are measured after compilation.
Data setInput →\to outputInput coord.nzn_{z}ddmmSpatial kernel
Forced isotropic turbulencef→u(⋅,0.5)\mathbf{f}\to\mathbf{u}(\cdot,0.5)21000337000Five-scale Matérn
Turbulent species transportinflow →u(⋅,0.5)\to\mathbf{u}(\cdot,0.5)1388237000Five-scale Matérn
Merging vorticesω(⋅,0)→u(⋅,0.4)\omega(\cdot,0)\to\mathbf{u}(\cdot,0.4)50032500Cubic PHS
Lid-driven cavityu(⋅,0)→u(⋅,5)\mathbf{u}(\cdot,0)\to\mathbf{u}(\cdot,5)2000121000Cubic PHS
Cylinder flow (no shedding)u(⋅,0)→u(⋅,10)\mathbf{u}(\cdot,0)\to\mathbf{u}(\cdot,10)2000121000Cubic PHS
Cylinder flow (vortex shedding)u(⋅,0)→u(⋅,10)\mathbf{u}(\cdot,0)\to\mathbf{u}(\cdot,10)2000121000Cubic PHS
Table 3: The ten benchmarks: input and output of the operator, number of sampled input coordinates, number nzn_{z} of principal input coordinates (the numerical rank of the training inputs, capped at three), spatial dimension dd , number mm of output points, and spatial kernel. Taylor–Green vortices (spacetime) has spacetime outputs, at the same 500500 points at four times.
Figure 6: All benchmarks at N=10000N=10000 . Left: test error of the fixed-shape kernel ( D∥=ID_{\parallel}=I , α⊥=1\alpha_{\perp}=1 ) without augmentation, divided by the test error of each other model (values above one are improvements). Right: density dKd_{K} of the selected kernel matrix.
Figure 7: Selected density dKd_{K} of the kernel matrix against the number NN of training functions for the three multivariate flows and the lid-driven cavity. The affine columns are not counted. Figure 21 shows the other benchmarks.
Figure 8: Test error (mean relative velocity error of the predictions for the 100 test inputs) against the number NN of training functions for the three multivariate flows and the lid-driven cavity. Dashed: with affine augmentation. Figure 20 shows the other benchmarks.
Data setNNOursPPKMVKMGeo-FNOTransolver
Forced isotropic turbulence100007.1×10−5\mathbf{7.1\times 10^{-5}}3.6×10−43.6\times 10^{-4}1.6×10−41.6\times 10^{-4}1.5×10−41.5\times 10^{-4}1.2×10−31.2\times 10^{-3}
Turbulent species transport100003.7×10−5\mathbf{3.7\times 10^{-5}}2.4×10−42.4\times 10^{-4}2.4×10−42.4\times 10^{-4}5.0×10−45.0\times 10^{-4}2.7×10−32.7\times 10^{-3}
Merging vortices5001.9×10−6\mathbf{1.9\times 10^{-6}}1.1×10−41.1\times 10^{-4}1.1×10−41.1\times 10^{-4}1.1×10−31.1\times 10^{-3}3.4×10−33.4\times 10^{-3}
Lid-driven cavity100004.7×10−7\mathbf{4.7\times 10^{-7}}1.4×10−61.4\times 10^{-6}1.5×10−61.5\times 10^{-6}1.5×10−41.5\times 10^{-4}2.1×10−42.1\times 10^{-4}
Cylinder flow (no shedding)1002.8×10−5\mathbf{2.8\times 10^{-5}}4.0×10−54.0\times 10^{-5}5.3×10−55.3\times 10^{-5}1.6×10−41.6\times 10^{-4}5.7×10−35.7\times 10^{-3}
Cylinder flow (vortex shedding)100003.0×10−7\mathbf{3.0\times 10^{-7}}3.1×10−73.1\times 10^{-7}2.3×10−62.3\times 10^{-6}3.7×10−53.7\times 10^{-5}1.5×10−41.5\times 10^{-4}
Table 4: Test error (mean relative velocity error of the predictions for held-out inputs) of our sparse local model (the one of the four with the smallest validation error) and of the methods reported by Sharma et al. [53] at the same training size NN : dense PPKM, the kernel method of Batlle et al. [2] on output samples (VKM), and Geo-FNO and Transolver, each with the better of its errors with and without a divergence penalty. The earlier errors are averages over all 200200 held-out pairs and ours over the 100100 test pairs. Our models of the cylinder flow without shedding at N=100N=100 were fitted for this comparison, with the settings of the other training sizes. Bold: smallest test error in each row.
Figure 9: Accounted storage of the fitted models (stored arrays and factors, not measured peak memory) against NN for the three multivariate flows and the lid-driven cavity. Figure 22 shows the other benchmarks.
Figure 10: Median time to predict all 100100 test fields after setup, including the kernel evaluation, the sparse solves, and the local output reconstruction. Fitting and setup are excluded (Appendix D.6.1 ).
Data setGlobal αspatial\alpha_{\rm spatial}Local αspatial\alpha_{\rm spatial}Largest gain (global)
Forced isotropic turbulence3.7×1093.7\times 10^{9}396\mathbf{396}1.21.2
Turbulent species transport2.9×1092.9\times 10^{9}9882\mathbf{9882}7.17.1
Merging vortices4.8×1094.8\times 10^{9}73\mathbf{73}868.1868.1
Lid-driven cavity5.8×1095.8\times 10^{9}13\mathbf{13}34.534.5
Cylinder flow (no shedding)5.1×1075.1\times 10^{7}149\mathbf{149}4.84.8
Cylinder flow (vortex shedding)5.2×1075.2\times 10^{7}149\mathbf{149}2.02.0
Table 5: Amplification factor αspatial=∥B∥F∥(A⋆)T∥F/∥B(A⋆)T∥F\alpha_{\rm spatial}=\|B\|_{F}\|(\mathbf{A}^{\star})^{\mathsf{T}}\|_{F}/\|B(\mathbf{A}^{\star})^{\mathsf{T}}\|_{F} of ( 66 ), which measures the cancellation among the spatial coefficients, for N=1000N=1000 training fields, for the global and for our local output reconstruction (stencils combined as in Section 5.3 ), and the largest ratio of the direct to the Lagrange error with global output reconstruction from Table 6 . Bold: smaller amplification factor in each row.
 Global output reconstructionLocal output reconstruction
Data setModelsMedianMaxModelsMax ∣| ratio −1∣-1|
Forced isotropic turbulence281.001.2285.9×10−115.9\times 10^{-11}
Turbulent species transport281.007.1288.6×10−78.6\times 10^{-7}
Merging vortices2848.96868.1284.2×10−84.2\times 10^{-8}
Lid-driven cavity285.0734.5281.4×10−71.4\times 10^{-7}
Cylinder flow (no shedding)281.054.8287.0×10−87.0\times 10^{-8}
Table 6: Ratio of the test error of direct evaluation to that of Lagrange evaluation, for the same unaugmented operator models (the fixed-shape kernel of Section 3.1 at the radii of the five prescribed densities and at its validated radius, and the fitted-shape kernel at its validated radius, for N=500N=500 , 10001000 , 50005000 , and 1000010000 ) with global and local output reconstruction. Values above one favor Lagrange evaluation. Over all 560 pairs, the smallest ratio is 0.999870.99987 .
Appendix figures & tables29 assets

Supplementary material from the paper’s appendix.

Appendix

Data setStencilsRetained rankTruncated stencils
Turbulent species transport700045–854338
Forced isotropic turbulence7000850
Buoyancy-driven cavity5000200
Lid-driven cavity1000200
Cylinder flow (no shedding)1000200
Cylinder flow (vortex shedding)1000200
Appendix
Table 7: Retained divergence-free polynomial ranks of the local systems on the nodes of each benchmark, at degree four and with periodic images, from the construction above. The full ranks are 20 in 2D and 85 in 3D. Taylor–Green vortices (spacetime) has the stencils of Taylor–Green vortices (space).
Figure 11: Second-order structure functions between the nodes of the reference fields and of the local PHS and multiscale Matérn output reconstructions of the test predictions at N=10000N=10000 .
Figure 12: Cylinder flow. (A) Initial speed. (B) Speed at t=10t=10 , on the same color scale.
Figure 13: Taylor–Green vortices. (A) Initial speed. (B) Speed at t=1t=1 , on the same color scale.
Figure 14: Lid-driven cavity. (A) Initial speed. (B) Speed at t=5t=5 , on the same color scale.
Figure 15: Backward-facing step. (A) Inflow profile. (B) Speed at t=5t=5 , with the transverse direction expanded.
Figure 16: Buoyancy-driven cavity. (A) Initial temperature. (B) Speed at t=5t=5 .
Figure 17: Merging vortices. (A) Initial vorticity. (B) Speed at t=0.4t=0.4 .
Figure 18: Turbulent species transport. (A) Inlet velocity on the plane z=0z=0 . (B) Speed at t=0.5t=0.5 on the whole domain.
Figure 19: Forced isotropic turbulence. (A) Forcing magnitude. (B) Speed at t=0.5t=0.5 .
NNDirectLagrangeDirect/LagrangeDirect on samples∥B∥F\|B\|_{F}
5001.54×10−51.54\times 10^{-5}1.20×10−5\mathbf{1.20\times 10^{-5}}1.281.20×10−5\mathbf{1.20\times 10^{-5}}1.58×1091.58\times 10^{9}
10008.39×10−68.39\times 10^{-6}4.25×10−64.25\times 10^{-6}1.974.24×10−6\mathbf{4.24\times 10^{-6}}2.24×1092.24\times 10^{9}
50003.19×10−63.19\times 10^{-6}9.54×10−79.54\times 10^{-7}3.349.32×10−7\mathbf{9.32\times 10^{-7}}5.00×1095.00\times 10^{9}
100003.42×10−63.42\times 10^{-6}6.64×10−76.64\times 10^{-7}5.146.42×10−7\mathbf{6.42\times 10^{-7}}7.08×1097.08\times 10^{9}
Appendix
Table 8: Test error at the nodes (mean relative velocity error over the 100 test inputs) of direct and Lagrange evaluation of the validated unaugmented fixed-shape model of the lid-driven cavity with global output reconstruction, of direct evaluation applied to the sampled velocities instead of the coefficients, and the Frobenius norm of the spatial coefficients of the training fields. Bold: smallest test error in each row.
Data setFixedFittedFixed + affineFitted + affine
Forced isotropic turbulence0.680.690.760.66
Turbulent species transport0.851.091.031.26
Merging vortices0.770.730.570.53
Lid-driven cavity0.971.090.961.08
Cylinder flow (no shedding)0.910.870.850.84
Cylinder flow (vortex shedding)0.540.530.440.47
Appendix
Table 9: Least-squares exponent τ\tau of the test error ∝N−τ\propto N^{-\tau} over N∈{500,1000,5000,10000}N\in\{500,1000,5000,10000\} . Bold: fastest decay in each row.
Data setFixedFittedFixed + affineFitted + affine
Forced isotropic turbulence2.80×10−32.80\times 10^{-3}1.62×10−31.62\times 10^{-3}9.77×10−59.77\times 10^{-5}7.14×10−5\mathbf{7.14\times 10^{-5}}
Turbulent species transport4.34×10−44.34\times 10^{-4}6.86×10−56.86\times 10^{-5}1.89×10−41.89\times 10^{-4}3.69×10−5\mathbf{3.69\times 10^{-5}}
Merging vortices3.02×10−7\mathbf{3.02\times 10^{-7}}3.07×10−73.07\times 10^{-7}3.39×10−73.39\times 10^{-7}3.19×10−73.19\times 10^{-7}
Lid-driven cavity6.42×10−76.42\times 10^{-7}4.76×10−74.76\times 10^{-7}6.42×10−76.42\times 10^{-7}4.72×10−7\mathbf{4.72\times 10^{-7}}
Cylinder flow (no shedding)2.02×10−72.02\times 10^{-7}1.94×10−71.94\times 10^{-7}1.70×10−7\mathbf{1.70\times 10^{-7}}1.73×10−71.73\times 10^{-7}
Cylinder flow (vortex shedding)3.19×10−73.19\times 10^{-7}2.95×10−7\mathbf{2.95\times 10^{-7}}3.19×10−73.19\times 10^{-7}2.96×10−72.96\times 10^{-7}
Appendix
Table 10: Test error (mean relative velocity error of the predictions for the 100 test inputs) at N=10000N=10000 . The two augmented models use the same affine columns. Each model is fitted and validated separately, and all four use the same local output reconstruction. Both Taylor–Green cases use the amplitude normalization ( 72 ). Bold: smallest test error in each row.
Data setFixedFittedFixed + affineFitted + affine
Forced isotropic turbulence10.680 / 235.210.680 / 234.010.685 / 244.110.682 / 233.4
Turbulent species transport2.662 / 126.63.130 / 126.02.630 / 123.33.137 / 137.4
Merging vortices1.397 / 79.61.029 / 62.40.772 / 37.70.767 / 39.0
Lid-driven cavity1.135 / 75.91.135 / 71.51.136 / 70.21.136 / 70.3
Cylinder flow (no shedding)1.141 / 73.81.141 / 73.71.142 / 75.71.142 / 74.0
Cylinder flow (vortex shedding)1.140 / 76.21.140 / 75.31.141 / 76.21.141 / 74.2
Appendix
Table 11: Accounted storage in decimal GB / median warm prediction time in ms for 100 fields at N=10000N=10000 . Storage includes the stored model arrays, the solve buffers, and cuDSS’s estimate of the memory of its factors. It is not measured peak GPU memory. Times include operator Lagrange evaluation and local output reconstruction. Bold: smallest storage and smallest time in each row.
Figure 20: Test error of the other six benchmarks. Dashed: with affine augmentation.
Figure 21: Selected density of the other six benchmarks. With augmentation, the Taylor–Green models select about 10%10\% density at the two largest NN , against 75%75\% or 90%90\% without it.
Figure 22: Accounted storage of the other six benchmarks (Appendix B.3 ).
Figure 23: Median warm time to predict all 100100 test fields for the other six benchmarks.
NNOperator modelTest errorDensity (%)Factor ( 10610^{6} )GBms
500Fixed shape0.021490.000.1250.46560.9
500Fitted shape0.01390.000.1250.46564.3
500Fixed shape + affine9.37×10−49.37\times 10^{-4}90.000.1250.46758.5
500Fitted shape + affine4.90×10−4\mathbf{4.90\times 10^{-4}}90.000.1260.46756.5
1000Fixed shape0.015790.000.9900.75064.8
1000Fitted shape0.013990.000.9910.75070.2
Appendix
Table 12: Forced isotropic turbulence: all models. Test error is the mean relative velocity error of the predictions for the 100 test inputs. Density is dK=nnz(KUU)/N2d_{K}=\mathrm{nnz}(K_{UU})/N^{2} in percent. Factor is the number of entries of the cuDSS Cholesky factor in millions. GB is the accounted storage and ms the median warm time for all 100 fields. Bold: smallest test error at each NN .
NNOperator modelTest errorDensity (%)Factor ( 10610^{6} )GBms
500Fixed shape4.39×10−34.39\times 10^{-3}63.940.1160.22654.3
500Fitted shape2.50×10−32.50\times 10^{-3}90.000.1290.22852.8
500Fixed shape + affine3.35×10−33.35\times 10^{-3}40.980.0970.22453.1
500Fitted shape + affine2.48×10−3\mathbf{2.48\times 10^{-3}}44.600.0960.22454.0
1000Fixed shape1.25×10−31.25\times 10^{-3}89.711.0000.27857.1
1000Fitted shape9.53×10−49.53\times 10^{-4}85.620.9620.27657.8
Appendix
Table 13: Turbulent species transport: all models. Test error is the mean relative velocity error of the predictions for the 100 test inputs. Density is dK=nnz(KUU)/N2d_{K}=\mathrm{nnz}(K_{UU})/N^{2} in percent. Factor is the number of entries of the cuDSS Cholesky factor in millions. GB is the accounted storage and ms the median warm time for all 100 fields. Bold: smallest test error at each NN .
NNOperator modelTest errorDensity (%)Factor ( 10610^{6} )GBms
500Fixed shape2.85×10−62.85\times 10^{-6}75.060.1040.0256.0
500Fitted shape2.51×10−62.51\times 10^{-6}75.000.1060.0256.8
500Fixed shape + affine1.94×10−61.94\times 10^{-6}37.360.0650.0217.9
500Fitted shape + affine1.66×10−6\mathbf{1.66\times 10^{-6}}33.640.0600.0216.9
1000Fixed shape1.41×10−61.41\times 10^{-6}63.810.3860.0607.3
1000Fitted shape1.45×10−61.45\times 10^{-6}60.310.3870.0588.6
Appendix
Table 14: Merging vortices: all models. Test error is the mean relative velocity error of the predictions for the 100 test inputs. Density is dK=nnz(KUU)/N2d_{K}=\mathrm{nnz}(K_{UU})/N^{2} in percent. Factor is the number of entries of the cuDSS Cholesky factor in millions. GB is the accounted storage and ms the median warm time for all 100 fields. Bold: smallest test error at each NN .
NNOperator modelTest errorDensity (%)Factor ( 10610^{6} )GBms
500Fixed shape1.20×10−51.20\times 10^{-5}22.850.0530.0469.2
500Fitted shape1.20×10−51.20\times 10^{-5}22.850.0530.04610.8
500Fixed shape + affine1.15×10−5\mathbf{1.15\times 10^{-5}}21.490.0480.04610.4
500Fitted shape + affine1.15×10−5\mathbf{1.15\times 10^{-5}}25.000.0550.0479.2
1000Fixed shape4.24×10−64.24\times 10^{-6}19.990.1900.09610.6
1000Fitted shape4.11×10−64.11\times 10^{-6}25.000.2230.09811.4
Appendix
Table 15: Lid-driven cavity: all models. Test error is the mean relative velocity error of the predictions for the 100 test inputs. Density is dK=nnz(KUU)/N2d_{K}=\mathrm{nnz}(K_{UU})/N^{2} in percent. Factor is the number of entries of the cuDSS Cholesky factor in millions. GB is the accounted storage and ms the median warm time for all 100 fields. Bold: smallest test error at each NN .
NNOperator modelTest errorDensity (%)Factor ( 10610^{6} )GBms
500Fixed shape2.91×10−62.91\times 10^{-6}89.110.1290.0539.8
500Fitted shape2.42×10−62.42\times 10^{-6}90.000.1280.0539.6
500Fixed shape + affine2.27×10−62.27\times 10^{-6}85.520.1280.05210.6
500Fitted shape + affine2.20×10−6\mathbf{2.20\times 10^{-6}}90.000.1280.05310.9
1000Fixed shape9.21×10−79.21\times 10^{-7}54.740.6740.11113.5
1000Fitted shape9.21×10−79.21\times 10^{-7}54.740.6740.11113.3
Appendix
Table 16: Cylinder flow (no shedding): all models. Test error is the mean relative velocity error of the predictions for the 100 test inputs. Density is dK=nnz(KUU)/N2d_{K}=\mathrm{nnz}(K_{UU})/N^{2} in percent. Factor is the number of entries of the cuDSS Cholesky factor in millions. GB is the accounted storage and ms the median warm time for all 100 fields. Bold: smallest test error at each NN .
NNOperator modelTest errorDensity (%)Factor ( 10610^{6} )GBms
500Fixed shape1.35×10−61.35\times 10^{-6}70.210.1220.05112.2
500Fitted shape1.35×10−61.35\times 10^{-6}70.210.1220.0519.7
500Fixed shape + affine1.11×10−61.11\times 10^{-6}68.610.1220.05112.6
500Fitted shape + affine1.09×10−6\mathbf{1.09\times 10^{-6}}75.000.1250.05112.0
1000Fixed shape1.19×10−61.19\times 10^{-6}27.920.2620.09911.6
1000Fitted shape8.54×10−78.54\times 10^{-7}25.000.2310.09811.9
Appendix
Table 17: Cylinder flow (vortex shedding): all models. Test error is the mean relative velocity error of the predictions for the 100 test inputs. Density is dK=nnz(KUU)/N2d_{K}=\mathrm{nnz}(K_{UU})/N^{2} in percent. Factor is the number of entries of the cuDSS Cholesky factor in millions. GB is the accounted storage and ms the median warm time for all 100 fields. Bold: smallest test error at each NN .
NNOperator modelTest errorDensity (%)Factor ( 10610^{6} )GBms
500Fixed shape9.41×10−79.41\times 10^{-7}46.170.0940.0286.5
500Fitted shape7.72×10−7\mathbf{7.72\times 10^{-7}}50.000.0970.0297.7
500Fixed shape + affine8.76×10−78.76\times 10^{-7}40.070.0830.0287.7
500Fitted shape + affine8.83×10−78.83\times 10^{-7}50.000.0970.0296.7
1000Fixed shape3.31×10−73.31\times 10^{-7}51.760.4220.0687.6
1000Fitted shape7.86×10−9\mathbf{7.86\times 10^{-9}}50.000.4140.0678.6
Appendix
Table 18: Taylor–Green vortices (space): all models. Test error is the mean relative velocity error of the predictions for the 100 test inputs. Density is dK=nnz(KUU)/N2d_{K}=\mathrm{nnz}(K_{UU})/N^{2} in percent. Factor is the number of entries of the cuDSS Cholesky factor in millions. GB is the accounted storage and ms the median warm time for all 100 fields. Bold: smallest test error at each NN .
NNOperator modelTest errorDensity (%)Factor ( 10610^{6} )GBms
500Fixed shape8.01×10−78.01\times 10^{-7}46.230.0930.0526.4
500Fitted shape6.58×10−7\mathbf{6.58\times 10^{-7}}50.000.0970.0536.3
500Fixed shape + affine7.38×10−77.38\times 10^{-7}39.510.0820.0528.1
500Fitted shape + affine7.65×10−77.65\times 10^{-7}50.000.0970.0536.8
1000Fixed shape2.71×10−72.71\times 10^{-7}51.580.4190.1167.9
1000Fitted shape9.30×10−99.30\times 10^{-9}50.000.4140.1159.8
Appendix
Table 19: Taylor–Green vortices (spacetime): all models. Test error is the mean relative velocity error of the predictions for the 100 test inputs. Density is dK=nnz(KUU)/N2d_{K}=\mathrm{nnz}(K_{UU})/N^{2} in percent. Factor is the number of entries of the cuDSS Cholesky factor in millions. GB is the accounted storage and ms the median warm time for all 100 fields. Bold: smallest test error at each NN .
NNOperator modelTest errorDensity (%)Factor ( 10610^{6} )GBms
500Fixed shape2.43×10−72.43\times 10^{-7}36.410.0730.0248.8
500Fitted shape2.43×10−72.43\times 10^{-7}36.410.0730.0248.6
500Fixed shape + affine9.46×10−8\mathbf{9.46\times 10^{-8}}34.940.0680.0249.6
500Fitted shape + affine9.46×10−8\mathbf{9.46\times 10^{-8}}34.940.0680.0248.6
1000Fixed shape9.31×10−89.31\times 10^{-8}49.590.3500.05911.6
1000Fitted shape1.14×10−71.14\times 10^{-7}25.000.2360.0508.4
Appendix
Table 20: Backward-facing step: all models. Test error is the mean relative velocity error of the predictions for the 100 test inputs. Density is dK=nnz(KUU)/N2d_{K}=\mathrm{nnz}(K_{UU})/N^{2} in percent. Factor is the number of entries of the cuDSS Cholesky factor in millions. GB is the accounted storage and ms the median warm time for all 100 fields. Bold: smallest test error at each NN .
NNOperator modelTest errorDensity (%)Factor ( 10610^{6} )GBms
500Fixed shape5.63×10−5\mathbf{5.63\times 10^{-5}}64.680.1180.15430.7
500Fitted shape5.63×10−5\mathbf{5.63\times 10^{-5}}64.680.1180.15429.6
500Fixed shape + affine8.18×10−58.18\times 10^{-5}10.000.0240.14829.1
500Fitted shape + affine7.39×10−57.39\times 10^{-5}10.000.0240.14828.6
1000Fixed shape3.03×10−5\mathbf{3.03\times 10^{-5}}17.330.1700.29832.5
1000Fitted shape3.03×10−5\mathbf{3.03\times 10^{-5}}17.330.1700.29832.2
Appendix
Table 21: Buoyancy-driven cavity: all models. Test error is the mean relative velocity error of the predictions for the 100 test inputs. Density is dK=nnz(KUU)/N2d_{K}=\mathrm{nnz}(K_{UU})/N^{2} in percent. Factor is the number of entries of the cuDSS Cholesky factor in millions. GB is the accounted storage and ms the median warm time for all 100 fields. Bold: smallest test error at each NN .
Data setOperator modelSelection (min)Rebuild (s)Initialization (s)
Forced isotropic turbulenceFixed shape3.8343.273.71
 Fitted shape26.3241.633.28
 Fixed shape + affine6.4543.704.15
 Fitted shape + affine27.0142.773.05
Turbulent species transportFixed shape3.241.484.17
 Fitted shape25.894.643.70
Appendix
Table 22: Offline and initialization costs at N=10000N=10000 . Selection includes kernel fitting, where applicable, and validation. Rebuild assembles and factors the selected operator. Initialization includes two synchronized 100-field prediction passes, compilation, and CUDA graph preparation. The four models share the setup of the output reconstruction and may reuse compiled code. The dash marks the selection time of the fixed-shape backward-facing-step model, which contains an external pause. Its raw time is in the released data. Bold: smallest time in each column for each data set.
More in this topic

PDE Operator Learning

Partial differential equation (PDE) operator learning builds models that map inputs such as coefficients and initial conditions to solution fields.

See all 228 papers

Explore similar work

CardsList
Feb 17, 2026

Fluids You Can Trust: Property-Preserving Operator Learning for Incompressible Flows

We present a novel property-preserving kernel-based operator learning method for incompressible flows governed by the incompressible Navier--Stokes equations. Traditional numerical solvers incur significant computational costs to respect incompressibility. Operator learning offers efficient surrogate models, but current neural operators fail to exactly enforce physical properties such as incompressibility, periodicity, and turbulence. Our kernel method maps input functions to expansion coefficients of output functions in a property-preserving kernel basis, ensuring that predicted velocity fields analytically and simultaneously preserve the aforementioned physical properties. We present universal approximation results and worst-case a priori convergence rates for our framework; empirically, the observed convergence rates exceed the pessimistic predictions across most benchmarks, motivating a formulation of more optimistic convergence rates. Another central contribution of this work is a novel computational framework revolving around streaming construction of kernel Gramians and recursive Schur-complement algorithms to solve the large block linear systems arising from property-preserving kernel methods. This kernel-based framework makes operator learning practical at scales up to 10,00010{,}000 training functions sampled at up to 10,00010{,}000 spatial locations. We evaluate the method on challenging 2D and 3D, laminar and turbulent, incompressible flow problems. Our method achieves up to six orders of magnitude lower relative ℓ2\ell_2 errors upon generalization and trains up to five orders of magnitude faster compared to neural operators. Moreover, while our method enforces incompressibility analytically, neural operators exhibit large divergence errors. Our results show that our method provides an accurate and efficient surrogate for incompressible flows.

Ramansh Sharma, Matthew Lowery, Houman Owhadi +1

Kahlert School of Computing University of Utah UT, USA · Department of Computing and Mathematical Sciences California Institute of Technology CA, USA

PDE Operator LearningOperator LearningComputational Fluid Dynamics
AskPDF
May 22, 2026

Operator Learning for Reconstructing Flow Fields from Sparse Measurements: a Language Model Approach

Reconstructing flow fields from sparse measurements is a fundamental problem in fluid mechanics with broad implications for modeling, control, and design. In this work, we propose a novel operator learning framework that leverages the architecture of language models to perform flow reconstruction in a mesh-free manner. We reformulate flow field reconstruction as a sequence-to-sequence learning task, where sparse measurements are treated as context and unobserved locations as queries. Our model learns to reconstruct the full flow field from sparse inputs, effectively capturing spatial correlations and long-range dependencies. We evaluate the proposed approach on four benchmark datasets: (1) two-dimensional vortex street simulations, (2) daily average temperature data across the contiguous United States, (3) three-dimensional blood flow simulations based on dissipative particle dynamics, and (4) three-dimensional turbulent jet flow measurements obtained via particle tracking velocimetry. Across all cases, our method demonstrates competitive reconstruction accuracy, even with highly incomplete data (less than 10% observed), and achieves efficient performance. The results highlight the potential of language models as robust and scalable tools for scientific data reconstruction, and suggest a promising direction toward the development of foundation models for scientific and engineering applications.

Qian Zhang, George Em Karniadakis

Division of Applied Mathematics, Brown University

PDE Operator LearningOperator LearningComputational Fluid Dynamics
AskPDF
Jul 7, 2026

Kernel-based Operator Learning: Error Analysis, Budget Allocation, and a Physics-Informed Extension

We study kernel-based operator learning in a two-stage sampling framework, where an offline kernel regression operator learns a discretized representation of the target operator from input-output pairs and an online kernel reconstruction operator recovers the output function from predicted observations. Our main theoretical contribution is an explicit budget allocation condition relating the number NN of training pairs, the number nn of input observations, and the output resolution mm. The condition is derived from a coupled error analysis that interprets the surrogate as a reconstruction from approximate data. This yields a decomposition of the total error into reconstruction and learning contributions that can be analyzed independently. As a consequence, we obtain quantitative scaling laws describing how NN, nn, and mm must be coupled to guarantee convergence and to balance offline learning and online reconstruction errors. The resulting estimates extend previous analyses of kernel-based operator learning. We further introduce a physics-informed extension that incorporates knowledge of the underlying PDE at evaluation time. Rather than encoding constraints directly into the kernel, we augment the online reconstruction step by penalizing PDE residuals at collocation points. The method requires no retraining for new inputs. Numerical experiments illustrate the theoretical findings and demonstrate the effectiveness of the proposed physics-informed reconstruction strategy.

Rüdiger Kempf

Applied and Numerical Analysis University of Bayreuth 95440 Bayreuth, Germany

Kernel MethodsPDE Operator LearningOperator Learning
AskPDF