computational engineers

Abstract: PhysicsNeMo-Mesh is the mesh module of the open-source PhysicsNeMo library. It provides mesh processing similar to PyVista, but on the GPU, in batches, and inside PyTorch's automatic differentiation. Like xarray for unstructured meshes, it attaches semantics to data: fields on points, cells, and whole meshes have names and shapes, and a simulation case keeps the complete PDE problem specification, including boundary conditions and their prescribed values. It addresses three problems that plague machine learning for physical simulation. First, GPU utilization is bottlenecked by I/O, because samples are very large relative to the compute per forward pass. Second, because traditional preprocessing is expensive, it is necessarily done once, during curation, rather than differently for each model and workflow at the point of use. Finally, advanced equivariant and mesh-native architectures require information that curation ordinarily discards. Fast GPU mesh processing makes it practical to preprocess at the point of use: on one GB200 GPU, mesh operations run $6.2\text{--}1{,}068\times$ faster than VTK with 72 threads, and loading a whole case is $8.0\text{--}88\times$ faster than reading its VTK files ($8.0\times$ when both read the same fields). Overlapped with training, this preprocessing costs $2.1\%$ of throughput for a $127\,\text{ms}$ Transolver step. So one dataset serves many models and workflows: in PhysicsNeMo's public recipe, one stored DrivAerML dataset feeds Transolver, GeoTransolver, and GLOBE. PhysicsNeMo-Mesh combines the speed of an intermediate array store with the flexibility of an archival format such as VTK. When shapes and wall conditions both vary, per-boundary conditions cut a heat-conduction error from $87\%$ to $6\text{--}9\%$. And by keeping field ranks and named boundaries, it retains the semantic structure that next-generation architectures will need.

Abstract: Neural-network-based PDE solvers require architectural choices, particularly the number of layers and neurons per layer. Fixing a single Gaussian radial basis function (RBF) feature layer and taking the number of features to the continuum limit removes the need to prescribe a finite feature count or select individual feature centers, thereby reducing representation arbitrariness. Motivated by this principle, we introduce Continuum-Feature Kernel Solver (CFKS), a deterministic framework for PDEs obtained by analytically marginalizing a bounded continuum of Gaussian RBF features. CFKS works directly with the resulting continuum-feature mean kernel and its operator-applied counterparts, parameterized by prescribed Gaussian widths and bounded center distributions. For linear PDEs, CFKS leads to symmetric operator-kernel collocation systems with null-space enforcement of physical constraints. Nonlinear PDEs are addressed through linearized initialization and nonlinear least-squares residual minimization using analytical kernel derivatives and Jacobians. Residual-guided collocation enrichment and sequential space--time-slab marching provide additional resolution strategies, while prescribed multiscale kernel refinement helps resolve localized nonlinear structures. Finite-feature experiments confirm convergence of empirical value and operator Gram matrices toward their continuum counterparts. Numerical studies include oscillatory and multiscale ODEs, irregular-domain Poisson and biharmonic problems, singular perturbations, transient transport equations, and viscous Burgers' equation. The results establish CFKS as a direct computational connection between finite Gaussian-RBF representations and continuum-feature kernel PDE solvers. Codes: https://github.com/vikas-dwivedi-2022/CFKS_Codes

Abstract: Learning physical models from spatiotemporal data requires simulating systems of partial differential equations (PDEs) and differentiating through their solution. Upon spatial discretization, these models give rise to stiff dynamical systems that require repeatedly solving large, sparse linear systems. Traditional scientific computing frameworks are largely CPU-based and require problem-specific tuning to accelerate, scale, and differentiate. In this work, we introduce diff-soup, a differentiable, sparse-oriented, unified PDE suite for GPU-native simulation and inference of large, stiff dynamical systems. Built in JAX, diff-soup provides 1) a fully differentiable design for problems with sparse structure, 2) scaling to large problems through the NVIDIA cuDSS sparse direct solver with model and data multi-GPU parallelism, and 3) end-to-end PDE-constrained optimization, from classical statistical inference to embedded deep learning. On stiff PDE benchmarks, diff-soup outperforms established CPU- and GPU-based frameworks by one to two orders of magnitude in both forward solves and gradients. We demonstrate both a classical Bayesian inference workflow, identifying constitutive relationships, and a modern deep learning workflow, training a neural field for the catalyst activity of a three-dimensional packed bed reactor. We provide diff-soup as an open-source repository to facilitate research in computational science and engineering.

Abstract: Walk-on-Spheres (WoS) enables mesh-free solutions of elliptic PDEs, but conventional methods repeatedly discard sampled trajectories, limiting their efficiency in interactive applications. We introduce reusable transport skeletons, a persistent representation that organizes stochastic walks into a shared computational graph. Our framework enables efficient reuse across spatial queries and boundary conditions, while supporting progressive refinement, structure-preserving compression, and local adaptation to dynamic geometry. It further enables interactive editing of transport between selected regions. Experiments demonstrate substantial computational savings while preserving solution accuracy. By transforming transient stochastic trajectories into reusable transport structures, our method provides a flexible foundation for interactive Monte Carlo PDE solving.

Abstract: The operation of intersecting two polyhedra, which can be reduced to clipping one polyhedron with a series of planes formed from another polyhedron, is central to intersection based mesh interpolation and material interface reconstruction. R3D is a robust library for such operation which implements the algorithm of Sugihara. In this paper, we present J3D, a new implementation of the clipping algorithm of Sugihara that provides substantial performance improvements over R3D.

Abstract: Reduced-order models simplify fluid simulation by representing complex flows with a small number of spatial modes. Building on Laplacian Eigenfluids, we introduce Neural Eigenfluid Bases (NEB), a neural representation whose basis functions can be evaluated anywhere in space. We enforce incompressibility and impermeable boundaries directly through differential lifts and potential boundary conditions. Geometry-local orthogonalization and Rayleigh minimization then learn an ordered basis. We connect learning in potential and velocity spaces, and derive a spectral-variance relation explaining why these objectives weight approximate modes differently, even though their exact positive spectra agree. Gauge fixing reduces the number of spatial derivatives needed for potential learning, while a Dirichlet energy loss learns independent circulation fields around obstacles. With 1024 modes, both objectives keep velocity errors below 1% in periodic vortex merger. For random flows around obstacles, potential-based learning captures more initial energy and gives lower final errors than numerical Laplacian bases of the same size. The reduced state is about a thousand times smaller than the full reference model. Three-dimensional tests separate spatial reconstruction errors from errors in coefficient evolution. Precomputed Galerkin operators allow millisecond-scale updates when the same basis is reused across simulations. This provides a physically constrained basis for continuous flow reconstruction and efficient reduced-order simulation.

Abstract: Hyperbolic partial differential equations (PDEs) transport information at finite speed along characteristics, often producing shocks and discontinuities that are challenging to capture numerically. Graph neural networks (GNNs) offer learned update rules, but neighbor aggregation alone does not tell a GNN how far information must travel and along which edge. As a result, when a GNN model takes a time step larger than what was used in training, shock fronts smear and oscillations grow. In this work, we ask whether making the scale of transport explicit as an edge-local Courant-Friedrichs-Lewy (CFL) signal can enhance GNNs' ability to model hyperbolic PDEs. Holding the architecture fixed, we vary only the information given to each edge, from a raw edge form to forms that progressively expose the edge-local transport interaction through direction, length, and wave speed. Because all of these are computed from the same ingredients, what changes is not the amount of information but its form. The results show a consistent pattern, in which feature count is not the limiting factor, and appending CFL directly to a raw edge form produces almost no benefit. The CFL conditioning proves most effective when the edge information already exposes the transport interaction, reducing final-step error by 10.4-14.4%, and the improvement widens as the step size grows larger. Notably, this improvement arises not from smearing the solution but from suppressing oscillations near shocks. These results point to a specific design principle that CFL conditioning is most effective when the edge representation already separates the direction of interaction from spatial scale, so that the CFL value has a well-defined quantity to modulate.

Abstract: PDE foundation models enable data-efficient adaptation but are deployed autoregressively, so short-horizon accuracy need not imply reliable simulation. We ask whether the short-horizon advantage of a foundation-model recipe persists when its predictions are repeatedly fed back as inputs for a long horizon. We define error amplification as G=E(H_f)/E(H_s), which measures full-horizon error relative to short-horizon error. We evaluate seven models on seven two-dimensional time-dependent PDE datasets: scratch-trained FNO and DeepONet; fine-tuned DPOT, GPhyT, Poseidon, and Walrus; and frozen in-context VICON. We report VRMSE, relative L_2, and relative MAE. Six datasets provide paired short- and full-horizon scores, while CE-RM contributes short-rollout sweeps only. At least one fine-tuned foundation model outperforms both scratch baselines in 26 of 30 short-rollout dataset--budget cells under every metric. Nevertheless, full-horizon VRMSE exceeds short-horizon VRMSE in 140 of 144 paired trained configurations, with a median error amplification of $2.04$. The count below VRMSE one falls from 104 to 53, and 12 full rollouts exceed ten although none does at the short horizon. Walrus achieves the lowest full-horizon VRMSE in 16 of 24 cells, yet its median error amplification is $2.70$. Relative L_2 and relative MAE show the same qualitative pattern. These results show that favorable short-horizon rankings can coexist with substantial error amplification. PDE foundation-model benchmarks should therefore report horizon-resolved absolute errors, amplification statistics, and failure incidence rather than a single local forecast score.

Abstract: Neuro-spectral architectures (NeuSA) are a class of physics-informed neural networks (PINNs) that have proven effective for solving time-dependent PDEs in the frequency domain. However, training and inference often rely on costly numerical integration of the learned vector field, typically using sequential Runge-Kutta time stepping, which limits temporal parallelism. We introduce PaPi, an integration scheme inspired by weakly nonlinear flows that enables time-Parallel inference through exponential integration and a Duhamel-Picard scheme. We show that the Duhamel-Picard iteration converges exponentially for all flows, and provide explicit error bounds that depend on how close the flow is to linear. Experiments show that PaPi accelerates training and inference by up to two orders of magnitude, achieving training times of under one second for low-dimensional PDE learning problems, while maintaining competitive predictive accuracy.

Abstract: Solving dynamical systems in complex and diverse geometries remains challenging for both numerical solvers and learning-based surrogate models. We propose Neural Operator Splitting, a framework for learning fluid dynamics while preserving the physical structure of fractional-step numerical schemes. Instead of learning the full-step evolution monolithically, we decompose each time step into physically interpretable steps, and learn each operator separately using intermediate states from the underlying numerical solver. Guided by the linearity properties of the fractional steps, we introduce a conditional-linear Fourier neural operator (CL-FNO) whose linear state mappings are modulated by nonlinear conditioning on inputs such as boundary geometry. We further introduce two complementary strategies, symmetry sampling and corrective cycling, to mitigate error accumulation during autoregressive rollout. Experiments on 2D and 3D fluid flows with varying geometries and initial conditions show that our full model reduces mean relative $L^2$ rollout error by 68\% in 2D and 40\% in 3D compared with a single full-step neural operator. At comparable error levels, neural-operator inference is approximately $10^3$ times faster per step than the Monte Carlo solver.

Abstract: Solutions of second-order elliptic PDE problems on nonconvex or cracked domains with mixed Dirichlet-Neumann conditions typically lie outside $H^2(Ω)$. Near re-entrant corners, crack tips and points where the boundary condition changes type, the gradient grows like $r^{λ-1}$ with $λ<1$. We study residual-based adaptive refinement in hybrid variational-collocation PINNs. A multiscale Fourier-feature network is trained on a composite loss. The loss combines a variational residual tested against Legendre bubble functions, a strong-form collocation term and penalties for the Dirichlet and Neumann data. The variational term needs only first derivatives of the network, whereas the collocation term requires second derivatives and provides a residual that can be evaluated at any point. After each training stage, candidate points selected by this residual are added to the collocation set. The method is tested on five two dimensional benchmarks, namely a Mode~III crack, a cross-shaped domain with four re-entrant corners with unequal amplitude, a reaction-diffusion problem with four localized peaks, a monotone semilinear problem with a boundary layer and the Motz problem. Compared with uniform point insertion under the same collocation and iteration budget, residual-based refinement reduced the relative $L^2$ error by up to a factor of $32$, the relative $L^\infty$ error by up to a factor of $38$ and the relative $H^1$ semi-norm error by up to a factor of $15$.

Abstract: We present KernelPack-JAX, a meshless, differentiable, GPU framework in JAX for solving forward and inverse PDE problems in domains whose boundaries, including internal obstacles, are given only by samples of their smooth pieces. Our node generator is based on a new stacked representation of piecewise-smooth boundaries, in which a parametric partition of unity of each smooth piece induces the supports of a curl-free implicit partition of unity of the whole boundary. This representation keeps the edges and corners of the boundary and classifies every point of space without a global fit. Boundary and interior nodes are then selected by a parallel greedy independent-set algorithm that reproduces sequential Poisson-disk sampling exactly. We discretize elliptic and parabolic problems with radial basis function-generated finite differences (RBF-FD) on stencils selected by a new visibility heuristic to avoid crossing the boundary, with ghost nodes for the boundary conditions. In addition, our framework exposes a GPU sparse direct solver to JAX as a differentiable linear solve, so that Jacobians of time-dependent problems cost a small multiple of a forward solve. We derive error estimates that separate consistency from stability and bound the parameter error by the forward error. We demonstrate high-order convergence for Darcy flow, reaction-diffusion, and advection-diffusion in 2D and 3D domains, including one with an internal obstacle, and show that the recovered parameters largely follow the forward errors. We also present timings for every stage.

Abstract: Multifidelity Monte Carlo (MFMC) reduces the cost of estimating the statistics of quantities of interest (QoIs) computed via expensive high-fidelity model by using cheaper low-fidelity models as control variates. We extend MFMC to vector-valued QoIs by replacing the scalar control-variate coefficients with matrices and minimizing a weighted mean-squared error. The optimal full-matrix coefficients are available in closed form, and their variance reduction is summarized by an effective correlation that takes the place of the squared Pearson correlation in the closed-form MFMC design. Restricting the coefficients to block-diagonal, diagonal, and scalar classes gives a nested hierarchy of estimators with increasing variance reduction and increasing pilot-data demand, with closed-form coefficients for compatible weightings. For a scalar QoI obtained by linearly projecting the state, correcting the full state before projecting is never worse at the population level. Experiments on a Gaussian problem, a reduced-order model of an advection--diffusion problem, and a coarse-grid two-phase flow problem show that the best class depends on the pilot size relative to the number of covariance quantities it must estimate, and that reusing the pilot samples lowers the mean-squared error.

Abstract: Residual-trained neural surrogates carry no error guarantee, a limitation for free-boundary problems such as phase change and contact. We derive guaranteed error bounds for neural approximations of parabolic obstacle problems in first-order state--flux--multiplier form. The main result bounds the flux error of the network and the energy error of its feasible projection $\max\{v,g\}$ by the constitutive residual, the initial mismatch, the complementarity defect, and the PDE residual weighted by $c_F^2/α$, without the coercivity and norm-equivalence constants of the first-order formulation. It extends the Apushkinskaya--Repin majorant for the heat operator to general coefficients, free multipliers and a flux bound. A penetration-consistency result shows that neural triples with arbitrarily small bound exist. Outward-rounded cellwise Taylor-model enclosures turn the bound into a finite computation valid on any adaptive partition. The flux and gradient errors of the projected output are certified below $0.27$ and $0.33$ of the solution scale for a one-dimensional Stefan problem and, after adaptive refinement, below one half for all three trained networks on a two-dimensional expanding disk. On an expanding ellipse the gradient bound stays above one half; substitution diagnostics trace the floor of the trained networks to the projected state near the free boundary, and the remaining gap to enclosure overhead. Without a closed-form solution, a validated lower bound replaces the solution norm; there the enclosure is nearly exact, but the trained networks have errors above the solution scale and no useful accuracy is certified. The bounds control a state--flux norm, not the multiplier or the free-boundary location.

Abstract: Physics-informed neural networks (PINNs) reduce or eliminate the need for large labelled datasets by enforcing governing equations during training. Physics-informed Deep Operator Networks (PI-DeepONets) extend this approach to learning solution operators, enabling rapid prediction across families of problems. Despite these advantages, we show that their application to higher-dimensional problems remains limited by the computational cost of training. To address this issue, we propose a modified PI-DeepONet architecture that reduces this cost by removing redundancy in how the Branch network processes input functions across collocation points. The proposed formulation is mathematically equivalent to the original architecture but reduces input tensor sizes and the number of Branch network evaluations. We analyse the computational complexity of both architectures and derive theoretical scaling laws for training time with respect to the number of collocation points and input functions. The proposed architecture is evaluated on 1D, 2D, and 3D linear advection problems. In all cases, it achieves comparable predictive accuracy to the original PI-DeepONet while consistently reducing training time. In 2D experiments, a maximum speed-up of $3.75\times$ is observed, with empirical training-time scaling agreeing with the predicted linear dependence on collocation points and input functions. Furthermore, training a 3D model on an NVIDIA A100 GPU is only feasible using the proposed architecture, as the original formulation exceeds memory limits. These results demonstrate that careful consideration of how PI-DeepONets are implemented can yield substantial gains, enhancing the practical performance of PI-DeepONets and enabling their application to otherwise intractable physics-informed learning problems.

Abstract: Randomized neural networks retain an output-linear least-squares solve, but globally sampled random features may fail to capture the spatial and scale structure of PDE solutions. We propose SFRBE--RaNN, a two-stage adaptive trial-space construction. SFR redistributes a fixed feature budget according to regional indicator masses and adapts regional weight scales through low-dimensional residual-based optimization. RBE constructs nonredundant enrichments by screening multiple activation families, removing represented response directions, and selecting the smallest admissible prefix. For solution-dependent operators, candidates are screened with a frozen operator and accepted only after nonlinear residual verification. For fixed linear problems, we establish an exact discrete residual-reduction identity and a conditional solution-error estimate. Building on this identity, we combine safeguarded SFR with RBE capture factors to derive a finite-horizon error bound involving baseline best approximation and empirical consistency, together with a conditional random-search success-probability bound. Numerical experiments on localized, multiscale, and nonlinear space--time problems demonstrate the method's effectiveness.

Abstract: This chapter provides a mathematical introduction to operator surrogates, which are finite-parametric maps approximating complex (typically nonlinear) operators between infinite-dimensional (function) spaces. A typical approach to constructing operator surrogates is the encoder-approximator-decoder paradigm, whereby infinite-dimensional inputs are mapped to a finite-dimensional latent space, processed, and then mapped back to the infinite-dimensional output space. We discuss various methods for creating efficient representation systems, including stable linear frameworks like Banach frames and data-driven techniques such as Principal Component Analysis (PCA), Active Subspaces, and nonlinear autoencoders. Subsequently, we survey key surrogate architectures, contrasting spectral methods, based on polynomial chaos expansions and related technologies, with more recent neural network approaches. We further discuss architectures that do not rely on representation systems, such as Fourier Neural Operators and transformers. We elaborate on their relationship and how they relate to other neural network operator surrogates, particularly those obtained by so-called "algorithm-unrolling". Additionally, we provide an introduction to the theoretical analysis of approximation and generalization errors. We explain how underlying low-dimensional structures, provided for example by holomorphy, are essential to overcome the so-called "curse of dimensionality": this refers to the fact that the complexity of numerical algorithms tends to increase exponentially in the problem dimensionality. We present convergence rate bounds for holomorphic operators that depend on the number of trainable parameters and the regularity of the input and output data. Throughout, we prioritize intuition over detailed technical proofs and provide precise pointers to the literature for detailed arguments.

Abstract: We present a method for enforcing exactly Dirichlet, Neumann, and Robin conditions on general triangular domains with arbitrarily curved boundaries. The algorithmic constructions developed here complements those in a previous work for general quadrilateral domains. Our method combines (i) the exact mapping between a general curved triangular domain and the standard triangle, (ii) TFC constrained expressions, and (iii) Barnhill transfinite interpolations, to provide explicit closed-form trial ansatze that mathematically satisfy the imposed boundary conditions exactly. Corner compatibility induced by the boundary conditions is pivotal, and enforcing exactly the imposed boundary conditions hinges on the exact enforcement of the induced compatibility constraints. The corner compatibility associated with the acute vertices of the standard triangle and the coupling between gradient components when the hypotenuse corresponds to a Neumann/Robin boundary pose particular challenges. Two strategies are critical to resolving these issues, by reformulating conditional corner constraints into equivalent unconditional forms and by introducing a separate TFC-style auxiliary edge construction together with an associated linear solve on the hypotenuse. We analyze the corner compatibility constraints and present formulations for the general form of functions that exactly satisfy the imposed Dirichlet, Neumann, or Robin conditions on curved triangular domains. The proposed method has been implemented using extreme learning machine type randomized neural networks for solving PDE problems. We present ample numerical experiments with linear/nonlinear PDE problems on a variety of domains with complex boundary geometries. Numerical results show that the current method has enforced the Dirichlet, Neumann, and Robin conditions on curved domain boundaries to the machine precision.

Abstract: This work presents a hybrid numerical approach for solving linear systems arising from the discretization of the two-dimensional parametric Helmholtz equation. A convolutional neural network based on the U-Net architecture is trained in an self-supervised manner (i.e., without explicitly generating reference solutions) to approximate the inverse of the discretized Helmholtz operator, using a loss function involving the residual norm of the linear system. The trained network is used as a nonlinear preconditioner within the flexible GMRES (FGMRES) algorithm. Numerical experiments show that while the neural network is not accurate enough to act as a standalone solver, it significantly improves the convergence of FGMRES when employed as a preconditioner. The neural preconditioner demonstrates robust performance and generalization capabilities with respect to variations in the velocity field and the domain size. Comparisons with algebraic preconditioners, based on LU factorizations of approximations of the original system matrices, indicate the superior efficiency of the neural approach under equivalent conditions. We believe that the proposed method is not tied to a specific neural architecture and can be extended to other parametric PDEs.

Abstract: To solve a family of large linear systems where the system matrix depends on the family parameter quadratically, we propose a new solution approach. Our method employs the derivatives of solutions with respect to the parameter. It turns out that the accumulated solution derivatives of increasing order, evaluated at a single specially chosen parameter value, provide successively better approximations to solution subspaces for the whole linear system family. We show that the derivatives can be computed by a recurrence, involving, at each step, solving a linear system with the same matrix. Constructing an orthonormal basis for the subspace spanned by the derivatives then leads to an efficient solver where just a single (sparse) LU factorization (or a single preconditioner setup in an iterative scheme) suffices to solve the whole linear system family. We call the new method ODB (orthonormal derivative basis method). As shown in numerical tests, compared to conventional solution methods, the ODB method allows to drastically reduce computational costs.

Abstract: Motivated by the success of large language models (LLMs) with broad generalizability and robustness to noisy or unreliable pre-training data, we explore in-context prediction of PDE solution fields using approximate pre-training data. In addition, inspired by the posterior predictive mean-based inference mechanism of the in-context learning in prior-data fitted networks (PFNs), we propose PDE-PFN, a prior-data fitted neural solver that approximates the posterior predictive mean of PDE solutions via in-context learning. PDE-PFN builds on a PFN architecture with self- and cross-attention mechanisms of Transformer and is pre-trained on noisy approximate solutions generated by physics-informed neural networks, serving as diverse but not necessarily exact priors. Through experiments on a range of two-dimensional PDEs, we demonstrate that PDE-PFN achieves empirical generalization across heterogeneous equations, robustness under noisy priors, and zero-gradient-update in-context inference capability. PDE-PFN achieves competitive or improved predictive accuracy relative to the evaluated baselines while supporting flexible context and query representations.

Abstract: Physics-informed machine learning has gained significant attention in recent years. In regimes of limited data, parametrized Gaussian processes have become popular. Existing approaches, however, often face limitations, such as requiring parametrizable (also called controllable) systems or a large number of output tasks. In this work, we introduce a systematic procedure we call "ghost tasking", using auxiliary tasks to circumvent these limitations. We prove that such ghost tasks can render any non-parametrizable system effectively parametrizable, enabling algorithmic construction of parametrized Gaussian Processes while keeping the number of required tasks (i.e. output dimensions) and latent functions low. We find that ghost tasking performs especially well in an inverse problem setting, even with very few available data. We show the usage and power of ghost tasking in three experiments, providing systematic comparisons to the only other currently available method applicable to all experiments. We provide necessary syntax and explications for two computer algebra programs that compute parametrizations for systems with polynomial or rational coefficients. Our theoretical results extend to systems with meromorphic functions.

Abstract: Physics-informed neural solvers often struggle with multiscale behavior and sharp gradients due to the use of global approximation spaces and nonconvex training. We introduce AB-PIELM, an adaptive-basis physics-informed extreme learning machine that combines domain decomposition with joint optimisation of spatial and spectral hyperparameters while retaining deterministic normal-equation training. The method adapts both the number and placement of subdomains and the locality of radial basis functions, enabling problem-dependent allocation of representational capacity without modifying the underlying solver. Experiments on oscillatory function approximation and singularly perturbed advection-diffusion equations demonstrate accurate solutions across a wide range of stiffness regimes while using substantially fewer neurons than existing neural and PIELM-based approaches. The learned hyperparameters are interpretable, concentrating resolution near sharp gradients. The framework also supports inverse problems, successfully recovering diffusion coefficients from sparse noisy observations using Bayesian optimisation. These results indicate that adaptive-basis PIELM formulations provide an efficient and stable alternative to gradient-trained neural PDE solvers and naturally extend to broader classes of differential equations.

Abstract: Assessing useful reuse in neural PDE solvers is challenging: final accuracy can reflect source learning and target-time computation. Our reuse contract separates solution accuracy, learning contribution, and numerical utility through paired state comparisons, matched target information and budgets, and cost accounting. A literature audit extracts 18 version-specific protocol records from 12 papers, documenting retained states, target-time resources, and reported controls. For a fixed linear system and residual tolerance, we construct two initial guesses with identical solution-error, energy-error, and residual norms, reaching the same solution with different conjugate-gradient (CG) iteration counts. Across 240 source-training trajectories, two linear PDE families, Fourier neural operators and convolutional networks, a fixed predictor's benefit reverses across correction algorithms. Among pairs with both relative prediction errors less than or equal to 5 percent on 64 in-distribution tasks (63 by 63 interior grids), reductions in all three norms accompany more CG iterations, at mean taskwise rates of 23.5 percent and 23.9 percent in two libraries. Work-based selection saves 2.50-3.33 CG iterations on held-out in-distribution tasks; matched adaptation demonstrates finite-budget pretraining value. Independent batches confirm a 0.73 percent complete online saving for one physics-trained Fourier neural operator against zero-initialized Poisson-preconditioned CG. Reuse requires matched state comparisons and downstream computational evidence.

Abstract: The Fourier Neural Operator (FNO) learns solution operators of partial differential equations (PDEs) through Fourier-space kernel parameterization, but frequency truncation can limit the learning of high-frequency variations. AM-FNO and SirenFNO generate kernels for all grid modes from spectral coordinates using shared networks, making coordinate encoding and generator design important. Recent work on implicit neural representations (INRs) has proposed constructing frequency interactions through explicit feature composition rather than relying on subsequent MLPs to form them implicitly. Building on this approach, we propose CAFE+FNO, which incorporates Content-Aware Frequency Encoding+ (CAFE+) into Fourier kernel generation. CAFE+ combines Fourier--Chebyshev features through parallel affine branches and a Hadamard product, forming interactions within and across the two feature families. A kernel MLP maps the resulting representation of each normalized spectral coordinate to a complex channel-mixing matrix. Each layer shares its generator across all stored modes, making the number of trainable parameters independent of the number of modes for a fixed architecture. We compare CAFE+FNO with existing FNO variants on five PDE benchmarks and conduct ablation studies on basis configuration, multiplicative composition, and bandwidth learnability. Code and experimental configurations are available at https://github.com/fabsk101/CAFEPlusFNO.git.

Abstract: We construct neural approximation classes native to the graph spaces H(div) and H(curl), in two and three dimensions and, for H(div), in any dimension. Every realization lies in the space for all parameter values, and with kinked potentials, such as ReLU networks, the admissible jumps appear at finite width. The classes are images of scalar and componentwise networks under fixed operators of the de Rham complex, and do not involve a mesh or finite element emulation. For H(div) in R^n two native classes are given on an equal footing, with a skew-symmetric potential $A$: $\mathrm{Div}\,A+R_nq+\mathbf{h}$, with the divergence $q$ as an explicit unknown, and $\mathrm{Div}\,A+\mathbf{z}$ with an $H^1$ field $\mathbf{z}$; for H(curl) the analogous classes are $\mathrm{grad}\,φ+Sr+\mathbf{h}$ in two dimensions and $\mathrm{grad}\,φ+\mathbf{z}$ in two and three dimensions. In all of them every interface jump of the field is carried by the potential term, $\mathrm{Div}\,A$ or $\mathrm{grad}\,φ$, while the remaining part has no interface jump (it is an $H^1$ field in the regular-decomposition classes); the classes with $\mathbf{z}$ are the componentwise approach enriched by this term. Known or learned interface geometry enters the potential through factors with trainable amplitudes, and the remaining part if the divergence jumps. The classes lead to the FOSLS-deRhaNN method, first-order system least squares with de Rham neural networks, whose loss is the least-squares functional posed in the natural spaces of the weak formulation; for elliptic equations this includes $H^{-1}$ right-hand sides and $H^{1/2}$ Dirichlet data. Elliptic equations with discontinuous coefficients and curl-curl problems are treated as instances, with the functional equivalent to the error; linear transport with discontinuous solutions and conservation laws with shocks use the same flux classes.

Abstract: Systems of equations arise in a wide range of scientific and engineering applications. The present work focuses on solvers for general systems of equations, including but not limited to those arising from partial differential equations. These systems can be broadly categorized into linear and nonlinear problems. For large linear systems, iterative solvers are generally preferred over direct methods due to the latter's superlinear growth of computational costs. Although convergence theory is well-developed under certain assumptions on the coefficient matrix, many classes of systems still pose open challenges. These difficulties become even more severe for systems of nonlinear equations, where nonlinear solvers typically rely on repeated linearization. For example, Newton's method may even converge quadratically near the solution; it can also converge slowly or diverge when the initial guess is not chosen appropriately. A wide range of solvers with diverse variants and hyperparameter settings exists, and the development of efficient and robust iterative methods remains an active area of research. Recently, machine learning (ML) techniques have been applied to enhance the efficiency of classical iterative methods while preserving their interpretability and reliability. We refer to these ML-enhanced iterative methods as hybrid iterative methods, in the sense that they combine classical iterative methods with ML. This paper provides a comprehensive overview of state-of-the-art approaches to constructing hybrid iterative methods for systems of both linear and nonlinear equations, while also discussing open challenges and outlining potential directions for future research.

Abstract: Neural solvers for partial differential equations (PDEs) can approach an accurate solution while their parameters grow without bound. In such cases, the limiting solution may have no finite representation in the chosen model, leaving the best loss unattained. Our analysis connects missing limits in deep neural tanh- networks to unbounded hidden parameters or increasingly redundant neurons. For a class of models built from translated kernels, we describe the missing functions and recover them by adding kernel derivatives to the model. This completion makes the best approximation attainable under standard assumptions. Numerical studies follow the associated parameter growth and explore how completion affects PDE optimization.

Abstract: Isogeometric analysis (IGA) solves partial differential equations accurately on exact NURBS geometry, whereas neural solvers are mesh-free but often orders of magnitude less accurate and typically trained by non-convex optimization without error control. We propose IGA-KAN, which uses local Kolmogorov-Arnold networks, fitted in closed form, to improve the IGA solution instead of replacing it. An IGA Galerkin solve produces u_h; on every knot-vertex patch a Kolmogorov-Arnold ridge model is fitted to the strong form of the equation, the exact boundary data and u_h, and the models are blended by IGA hat functions. With fixed inner functions the fit is one batched linear least-squares problem, without optimizer, learning rate or initialization. An a posteriori safeguard, motivated by a maximum-principle bound, decides where local models are used, keeping the IGA solution elsewhere. On eight benchmarks with exact solutions, five from the literature and one also posed on a domain fitted to a brain slice from MRI, the method reduces the error of IGA, at an unchanged number of Galerkin unknowns, by factors of 4.2 to 90 in L^2 and 4.1 to 220 in H^1 on the reference meshes, and its L^2 error is 6 to 6x10^4 times smaller than that of the best Kolmogorov-Arnold network trained from scratch on the same equations with a fixed budget. In an inverse problem it recovers an unknown constant source from one noise-free observation 167 times more accurately than IGA. The gain is attributed to the superconvergence of local averages of the Galerkin solution.

Abstract: Graph Neural Operators (GNOs) provide flexible surrogate models for learning solution operators of partial differential equations (PDEs). However, standard GNOs typically parameterize the integral kernel with a monolithic neural network and evaluate kernel interactions over graph edges, leading to substantial computational and memory overhead at high resolutions or with large neighborhoods. To address these limitations, we propose Functional Tensor Decomposition Graph Neural Operator (FTD-GNO), a memory-efficient GNO framework that decouples the high-dimensional continuous integral kernel into low-dimensional mode-wise functions. By instantiating the kernel with classical tensor decomposition formats, including CP, Tensor-Train, and Tucker decompositions, FTD-GNO enables algebraic reconstruction of the integral operator without explicitly materializing full edge-wise kernel tensors. This factorized formulation reduces the memory footprint of kernel evaluation and aggregation while retaining the continuous operator-learning structure of GNOs. Theoretical complexity analysis shows that FTD-GNO substantially lowers parameter and activation-memory costs associated with high-dimensional kernel construction. Experiments show lower peak memory than the corresponding unfactorized graph-integral baselines, with shorter recorded training times. Fourier-graph experiments further demonstrate that FTD can improve the efficiency of a graph-integral layer within a hybrid operator and has good scalability.

Abstract: The quadratic attention cost of Vision Transformers (ViTs) forces a trade-off between spatial resolution and rollout horizon, particularly for fine-scale PDEs where shocks, reaction fronts, and material interfaces occupy small, evolving regions of the domain. Conventional neural surrogates also lack mechanisms to adapt resolution dynamically. We propose WAMRViT, a ViT that tokenizes inputs as balanced quadtrees using a wavelet-inspired refinement criterion, jointly encodes position and refinement level with 3D rotary positional embeddings, and regrids in cell space during inference for stable long-horizon rollouts. A multi-scale variant retains each leaf at its native source resolution and lets the model learn across resolution levels. Unlike fixed-budget adaptive-tokenization methods, WAMRViT imposes no predetermined token count and supports fully adaptive topology throughout autoregressive rollout. To our knowledge, it is the first machine-learning surrogate to natively tokenize multi-level Adaptive Mesh Refinement (AMR) data. On uniform-grid benchmarks, uniform-patch WAMRViT improves finest-level region-of-interest VRMSE over a finest-patch uniform ViT while using substantially fewer tokens. The multi-scale variant achieves the lowest first-step full-field RMSE and VRMSE on both benchmarks and improves rollout-averaged full-field and refined-region accuracy at long horizons. With parallelized regridding, its end-to-end rollout cost lies between finest-patch and approximately token-matched coarser-patch ViTs. On a complex AMR combustion problem whose finest features cannot be represented natively by the evaluated uniform-grid baselines, WAMRViT operates directly on adaptive cells and substantially reduces finest-level error at matched transformer capacity. Code: https://github.com/tonyzyl/wamrvit

Abstract: Learning partial differential equation (PDE) dynamics across varying domains is central to predictive modelling and data-driven discovery of governing equations. However, geometric variation alters both field representation and the governing differential operators, confounding geometric effects with intrinsic physical properties in the observed dynamics. This work identifies geometry-physics confounding as a unified failure mechanism for PDE learning across varying domains. In forward operator learning, this confounding increases the burden of inferring geometry-dependent operator changes from finite data, reducing data efficiency and generalisation. In equation discovery, omitting geometry-induced operators misspecifies the candidate library, leading to biased parameters, missed governing terms and spurious terms. We propose a de-confounding framework that makes the known geometry-to-operator transformation explicit. Geometry-induced coefficient fields improve prediction and data efficiency across five operator-learning benchmarks, while geometry-complete candidate libraries recover the generating equations and reduce held-out PDE residuals by more than two orders of magnitude in both evolving-domain systems. By separating known geometric action from intrinsic physics, the proposed framework supports more reliable and data-efficient PDE learning across scientific and engineering problems with varying geometries.

Abstract: Neural solvers offer efficient surrogates for numerical simulation of partial differential equations (PDEs). For time-dependent problems, strong one-step accuracy does not necessarily translate into reliable autoregressive rollout. We observe that a solver based only on physical-state modeling can achieve lower one-step error, whereas its spectral-only counterpart can become more accurate at later rollout steps. Motivated by this observation, we present Transolver-$σ$, a neural PDE solver based on joint spectral--physical subspace modeling. Within each block, adaptive physical-state interactions and spectral transformations are modeled in dedicated latent subspaces, whose responses are recomposed to enable information exchange between the two representations. Within the physical subspace, we introduce Slice-Residual Physics-Attention (SRPA), which preserves an explicit slice-space identity path while retaining learnable cross-slice interaction. In parallel, an axis-factorized Fourier operator captures global spectral structure. Across five well-established PDE benchmarks spanning steady-state prediction and time-dependent dynamics, Transolver-$σ$ achieves state-of-the-art with a benchmark-averaged relative error reduction of 33.4% over the strongest baseline for each metric, while consistently improving autoregressive rollout over single-operator counterparts. Transolver-$σ$ further delivers strong gains on coupled multiphysics systems and real-world fluid and combustion measurements from RealPDEBench, demonstrating its effectiveness beyond standard simulation benchmarks.

Abstract: Physics-informed neural networks (PINNs) solve partial differential equations (PDEs) by incorporating governing physical laws into the training loss. For evolution equations, however, their conventional pointwise space--time representation does not explicitly encode temporal dependence, which can hinder accurate prediction. To mitigate this limitation, this work proposes a novel neural architecture termed a causal integral neural network (CinNet). The core module of CinNet is a Volterra-type causal integral term, which aggregates historical features to encode temporal dependence, thereby incorporating temporal causality at the architectural level rather than through training-level modifications as in many existing methods. Building on CinNet, we further develop a causal integral physics-informed neural network (CI-PINN) for solving evolution equations. Extensive numerical experiments on benchmark evolution equations demonstrate that the presented method outperforms various baseline PINN variants in terms of solution accuracy, with pronounced superiority under sparse-collocation scenarios. Additional empirical analyses show that CI-PINN exhibits low sensitivity to hyperparameter choices, while ablation studies confirm the effectiveness of the proposed network components.

Abstract: Neural operators are typically trained in a supervised fashion, which requires a dataset to be generated with a classical solver. Training them physics-informed, i.e., purely from the governing equations, removes this large offline cost and allows fresh samples to be drawn at every optimization step, but has so far been limited to simplified problems and trails supervised training in accuracy. The obstacle is the ill-conditioning of physics-informed losses, which differential operators induce and which worsens as the discretization is refined. We therefore propose a preconditioned residual loss function and show mesh-independent conditioning for elliptic problems and greatly improved conditioning for saddle point problems. Realized through geometric and algebraic multigrid, the construction applies to linear and nonlinear equations, steady or time-dependent, on structured and unstructured meshes, is agnostic to the neural operator architecture, and adds no cost at inference. On the Poisson, Allen-Cahn and stationary Stokes equations, the resulting label-free training matches supervised training and is four to twenty-five times more accurate than previous physics-informed operator learning methods.

Abstract: Accurate numerical solutions of partial differential equations (PDEs) are crucial in numerous science and engineering applications. In this work, we introduce a novel neural PDE solver named AFDONet, which incorporates neural operator learning and adaptive Fourier decomposition (AFD) theory for the first time into a specifically designed variational autoencoder (VAE) structure, to solve a general class of nonlinear PDEs on smooth manifolds. AFDONet is the first neural PDE solver whose architectural and component design is fully guided by an established mathematical framework (in this case, AFD theory), turning neural operator design from an art to a science. Thus, AFDONet also exhibits exceptional mathematical explainability and groundness, and enjoys several desired properties. Furthermore, AFDONet achieves outstanding solution accuracy and competitive computational efficiency in several benchmark problems. In particular, thanks to its deep connections with AFD theory, AFDONet shows superior performance in solving PDEs on i) arbitrary (Riemannian) manifolds, and ii) datasets with sharp gradients. Overall, this work presents a new paradigm for designing explainable neural operator frameworks.

Abstract: Automatic differentiation (AD) lets neural networks compute derivatives of governing equations to machine precision, and this precision has made it the computational backbone of physics-informed machine learning. Yet exactness in the mathematical sense is not the same as fidelity to the physics. Here I argue that a derivative can be numerically perfect and still be the wrong derivative for the problem at hand, because AD, by construction, has no notion of the physical structure a solution must obey. Convection and its associated directionality, diffusion, and dispersion are only the most visible instances of a much longer list that spans all branches of computational science and engineering, including conservation, thermodynamic consistency, symmetry, symplectic structure, positivity, monotonicity, and boundedness. Recognizing this broader gap reframes how the field should build the next generation of PDE-driven neural surrogates.

Abstract: Prior information for real-world physical quantities is most elegantly expressed via partial differential equations (PDEs). In this paper, we propose a novel way to solve PDEs using probabilistic inference on a factor graph. In general, factor graphs provide a natural way to encode prior knowledge into a model as explicit factors; here, this knowledge is provided by a governing PDE, which narrows the solution space, while observed data further shape the posterior over the parameters. The approximate parameter posterior is inferred using message passing based on moment matching, without posterior sampling or global gradient-based optimization. We demonstrate our approach on the first-order advection and the second-order semi-linear Fisher-KPP equations, where it achieves predictive accuracy comparable to a standard baseline while providing structured predictive uncertainty. Moreover, the inferred posterior marginal means and uncertainty structure match more closely those obtained using Hamiltonian Monte Carlo than the evaluated mean-field variational inference baseline, while requiring up to 10x less training time in our experiments, with inference speed comparable to variational inference.

Abstract: Physics-informed neural networks (PINNs) provide a learning-based framework for solving partial differential equations (PDEs), yet their training behavior can change substantially throughout optimization. Residual distributions, gradient interactions, regional learning difficulty, and model-capacity requirements may evolve over time, while the network architecture and major training mechanisms are typically determined before training. We propose PINNMorph, an online PINN adaptation framework based on large language model (LLM)-guided policy evolution. PINNMorph maintains a population of state-conditioned adaptation policies that map execution diagnostics to controlled interventions over topology modification, additive representation augmentation, objective balancing, gradient handling, adaptive sampling, and optimizer-phase control. At each intervention opportunity, candidate programs are instantiated from the current policy population, selected according to the observed training state, and applied directly to the PINN under training. The resulting model inherits its existing parameters and training state and continues optimization along the same trajectory. Execution outcomes are subsequently used to evaluate interventions and evolve the policy population. Unlike pre-training architecture search or fixed adaptation rules, PINNMorph jointly adapts the current PINN and the policies governing its interventions using feedback from actual training. Experiments on 13 PDE benchmarks show that PINNMorph achieves lower solution errors than SA-PINN, ConFIG, RoPINN, HARMONIC, and PINNsAgent across all evaluated problems. Ablation studies further examine the effects of online adaptation, state-conditioned intervention selection, and execution-feedback-driven policy evolution.

Abstract: We formalise in Lean 4, on top of Mathlib, the solvability of the Dirichlet problem for second-order linear elliptic operators in divergence form. The machine-verified results, with no sorry in the development, include the Poincaré inequality, the existence of weak solutions by the Lax-Milgram theorem, Rellich-Kondrachov compactness, the Fredholm alternative, the spectral theorem, interior regularity estimates, and the Sobolev embedding theorem. From these results we obtain a formalisation of classical solvability for sufficiently regular coefficients and data. Our Lean library includes a self-contained theory of Sobolev spaces developed independently of existing formalisations. Throughout the paper we associate each prose statement with the named machine-checked Lean declaration that discharges it.

Abstract: We investigate the use of neural network solvers for infinity and $p$-Laplace problems, which are fundamental in nonlinear analysis and have practical applications. Our approach employs Physics-Informed Neural Networks (PINNs) and Deep Operator Networks (DeepONets) to address computational challenges associated with large $p$ values, ranging from $2$ to $1000$, on various 2D and 3D domains. Our method offers advantages over traditional physics-based solvers, especially in three dimensions where mesh-based solvers become very costly for these problems. We also establish conditional convergence results for PINN approximations of both problems and a universal approximation result for DeepONet on the parametric $p$-Poisson problem. We demonstrate the effectiveness of these neural network solvers through numerical experiments and compare their performance with conventional methods.

Abstract: We propose the Memory-Efficient Neural Operator (MENO) as a high-performance PDE neural solver based on the Manifold Function Encoder (MFE). MENO features three primary advantages: (1) MENO has a significantly smaller memory footprint and much faster training speed than other popular architectures, with the memory footprint being independent of the data resolution, and therefore holds the potential for scaling up to large-scale models. (2) MENO can accept PDE inputs of arbitrary form, including arbitrary geometric domains and arbitrary discretizations. In particular, it is capable of handling cross-geometry scenarios, i.e., where the input functions and the output solutions are defined on different manifolds. (3) MENO exhibits strong generalization capability, and achieves the best accuracy on most of the benchmarks we tested, compared with the results reported in the literature. The code is available on GitHub at https://github.com/jpzxshi/MENO, and all numerical examples in this paper can be run with a single command to reproduce the reported results.

Abstract: In this work, we propose a hybrid iterative deep Ritz method (H-IDRM) for a class of interface problems for second-order elliptic operators. It is based on a new mixed formulation of the problem and involves solving a sequence of convex minimization problems. We employ a level-set neural network architecture, featuring a level-set representation of the interface, to accommodate the piecewise smoothness of the solution and the flux. The approach involves only volumetric representations instead of duality pairing on the interface and avoids explicit interface sampling that is inconvenient for complex interface geometries. Further, we present an analysis of the method, including the errors arising from the neural network approximation, Monte Carlo approximation, iterative scheme, and penalty parameters. Numerical experiments indicate that the H-IDRM outperforms existing neural solvers on problems with high-dimensional domains, intricate interface geometries, and lower subdomain regularity.

Abstract: Green's function based operator approximation offers an efficient route for solving linear partial differential equations under varying boundary conditions and source terms. Once the Green's function is learned, solutions for new configurations are obtained through integration rather than by solving the differential equation again. Existing Green's function learning methods face two structural challenges. The first is the singular behavior of the Green's function near the source point, which places a difficult approximation burden on neural networks. The second is the absence of explicit consistency between the learned Green's function and its gradient, although both quantities enter the integral solution representation directly. This work proposes SDC-GON, a Singular Decomposition and Consistency-Regularized Green's Operator Network that addresses both challenges within a unified framework. The Green's function is decomposed into an analytically known singular component and a smooth correction learned by the network, so that the neural approximation targets only the regular part of the response kernel. A self-consistency loss enforces agreement between the gradient and the autodifferentiation gradient of the smooth correction. The method is evaluated on two dimensional Poisson, three dimensional heat conduction, heterogeneous reaction diffusion, and Stokes benchmarks, consistently outperforming the compared baselines across all cases. On the heterogeneous pipe benchmark, SDC-GON achieves a testing error of $3.70\times10^{-4}$ with a smaller network architecture, compared with $9.60\times10^{-4}$ for the same-width baseline and $4.63\times10^{-4}$ for a larger configuration, demonstrating that structural improvements are more effective than increasing model size.

Abstract: Physics-informed neural networks (PINNs) require coordinated choices over network representation, sampling, loss construction, and optimization, while effective configurations often vary substantially across partial differential equations (PDEs). Existing automated PINN design methods can search candidate configurations, but information revealed during actual training is still used mainly for evaluation rather than to improve subsequent design, leading to repeated trial-and-error and inefficient use of training budget. We propose PINNsForge, an LLM-driven evolutionary framework for execution-feedback-based automated PINN design. PINNsForge generates diverse candidate configurations from PDE-related prior knowledge, evaluates them through actual training, and feeds high-performing designs together with accumulated execution evidence back to the LLM. Guided by observed optimization behavior, the LLM then refines, recombines, and explores coupled PINN design components, forming a continual cycle of generation, execution, feedback, and evolution. Unlike one-shot search or evaluation-only feedback, PINNsForge progressively converts training experience into improved design decisions for the target PDE. Across 25 PDE benchmarks, PINNsForge achieves the lowest mean MSE on 24 tasks compared with RoPINN, PINNsFormer, and PINNsAgent. Ablation studies further confirm the importance of the PDE knowledge base, execution feedback, and evolutionary search: removing these components increases the mean MSE to 3.74$\times$, 12.10$\times$, and 10.10$\times$ that of the full PINNsForge, respectively.

Abstract: Amortizing physics-informed neural networks (PINNs) across related PDEs requires describing each equation to a reusable solver. Coefficient vectors encode numerical parameters in predefined slots, leaving operator and cross-field assignments implicit. We make these relationships explicit in an operator graph, with nodes for fields, derivatives, terms, and residuals and coefficients retained as term attributes. A graph hypernetwork generates diagonal codes that initialize a meta-trained factorized PINN for each target equation. Meta-training and target-specific adaptation use governing equations and prescribed conditions without solution labels. We compare coefficient-vector, DeepSets-based term-set, and graph conditioning by solution accuracy within a fixed adaptation budget. In scalar convection-diffusion-reaction problems, both term-based descriptors improve high-reaction accuracy, with similar performance. In two-field Fisher-KPP, meta-training sees uncoupled and one-way systems; after 3,000 adaptation steps on unseen two-way coupling, the graph's mean final error is 35.7% below the term set and 67.7% below the coefficient vector. In a fixed-structure capacitively coupled plasma model, the coefficient vector performs best. These results support extending coefficient conditioning with explicit equation relationships for physics-based solver adaptation.

Abstract: Hybrid domain decomposition methods provide a flexible framework for coupling full order models (FOMs) and reduced order models (ROMs), but typically assume the model assigned to each subdomain is fixed throughout a simulation. This is limiting for transient problems in which localized features propagate through the domain and the regions requiring high-fidelity resolution change over time. We introduce a reinforcement learning (RL)-based approach for online adaptation of FOM-ROM models coupled via the overlapping Schwarz alternating method (O-SAM), an iterative domain decomposition method that solves subdomain-local problems while exchanging solution information through transmission boundary conditions on overlapping interfaces. Deep Q-networks (DQNs) are trained offline to select among subdomain-local FOMs and pre-trained Operator Inference (OpInf) ROMs using a reward balancing accuracy, cost, and model-switching frequency. Once trained, the policies are deployed predictively on problem instances not seen during training, without requiring a reference FOM solution. We demonstrate the approach on two examples: a 1D advection-diffusion problem with a moving front, and a 3D linear elastic wave propagation problem implemented in the Norma.jl solid mechanics code. For the advection-diffusion benchmark, the learned policy dynamically allocates high-fidelity resolution as the front propagates and outperforms static FOM/ROM assignments; letting the agent also adapt the domain decomposition provides no further benefit. For the elastic wave benchmark, learned policies for two and three subdomain decompositions track the propagating wave by assigning FOMs to subdomains containing the wave and ROMs elsewhere, as expected. Our results demonstrate the potential of RL to enable predictive online adaptation of model fidelity within Schwarz-based hybrid simulations.

Abstract: We introduce a novel online adaptive non-intrusive reduced-order modeling strategy for parameterized dynamical systems involving parameter and time-dependent reduced bases. The proposed framework is based on a unified Grassmann manifold formulation combining three key components: interpolation of local reduced subspaces for unseen parameters, geodesic online subspace updates driven by incoming high-fidelity snapshots, and a latent-space regression strategy relying on Grassmann-distance weighting and Procrustes alignment to consistently aggregate predictions from multiple local models. The adaptive reduced-order model is embedded in a partitioned fluid-structure interaction framework, where it predicts fluid interface forces to provide accurate initial guesses for the nonlinear coupling iterations, thus achieving computational speedups with no loss of accuracy. The reduced basis and the regression operators are adapted independently during the simulation and without requiring the storage of high-dimensional streaming data, preserving computational efficiency while substantially improving predictive capabilities. Numerical results on reference FSI test cases demonstrate superior accuracy with respect to static and global reduced-order models, leading to a significant reduction in the number of fixed-point iterations required for convergence. The proposed framework offers a flexible and fully non-intrusive approach for the efficient simulation of nonlinear parameter-dependent multiphysics problems.

Abstract: Operator learning for parametric partial differential equations (PDEs) typically builds global models over prescribed domains, requiring cross-condition data or costly physics-constrained training. Here we show that the output Jacobian of a neural solution model trained at one condition defines a reusable response space for cross-condition solution variations. We introduce Linearized Subspace Transfer (LST) to exploit this space and recover target solutions by minimizing the target PDE-system residual over response-space coordinates. Because any single response space has finite coverage, Active Transfer Modeling (ATM) uses post-transfer residuals as coverage indicators to selectively acquire response spaces from additional single-condition models. Across six systems, single-condition response spaces supported cross-condition transfer, with enrichment improving accuracy when added spaces expanded representation capacity. Relative to evaluated physics-informed operator baselines, ATM reduced error and offline construction cost, with orders-of-magnitude accuracy gains in representative cases and millisecond-to-second target adaptation. These results establish neural solvers as reusable local parametric models.

Abstract: We propose a linearized Physics-Informed Neural Network (lPINN), a reduced-order neural basis method for forward and inverse differential equations. In an offline stage, lPINN learns operator-compatible continuous neural basis functions from an ensemble of numerical solutions. The basis functions are differentiable through automatic differentiation and are pretrained using solution data together with either derivative information or physics residuals. For each new problem instance, the basis functions are frozen and the solution is obtained by minimizing the governing-equation residual together with applicable initial, boundary, regularization, and observational terms. Unlike surrogate and operator-learning methods, the training data define the trial space offline, while the instance-specific solution is computed online by enforcing the governing physics. Relative to vanilla PINNs, lPINN pretrains the nonlinear hidden-layer representation offline and performs online inference only in the final linear layer. We evaluate lPINN on forward and inverse problems for the advection-diffusion equation, Burgers' equation, and the nonlinear pendulum equation. Compared with vanilla PINNs, lPINN achieves lower solution and parameter errors while reducing online inference times by approximately one to more than three orders of magnitude, with the largest gains generally observed for limited residual or measurement data. Cross-resolution experiments show that the learned continuous representation can be evaluated on finer meshes without retraining and with nearly unchanged accuracy.

Abstract: Scientific machine learning methods such as physics-informed neural networks (PINNs) increasingly rely on domain decomposition for better scalability while solving partial differential equations (PDEs) over complex geometries, yet the resulting composite loss comprising residual, boundary, and interface terms is highly susceptible to conflicting gradients that degrade training. This work bridges domain decomposition with projection-based gradient surgery to systematically mitigate such conflicts in 2D and 3D settings. We evaluate two existing projection-based algorithms, PCGrad and ConFIG, and identify their performance degradation in specific scenarios such as 3D domains with multiple overlapping interfaces. To address this limitation, we propose Norm-PCGrad, a normalized variant that achieves state-of-the-art accuracy across a range of 2D and 3D domain decomposition problems. Across the benchmarks considered, Norm-PCGrad consistently achieves the lowest relative $L_2$ error compared to training without gradient surgery as well as to existing algorithms such as PCGrad and ConFIG, while incurring negligible additional computational overhead. To improve computational efficiency of domain decomposition frameworks such as Extended PINN (XPINN), we propose replacing vanilla PINNs in selected subdomains with separable architectures such as Separable PINN (SPINN), reducing the computational cost from quadratic (or cubic) to linear. We additionally demonstrate that gradient surgery extends to physics-informed Kolmogorov-Arnold Networks (PIKANs), yielding substantial accuracy improvements for 3D domain decomposition and confirming the generality of the proposed approach across network architectures.

Abstract: Parametric PDEs can admit different boundary conditions with different accuracy and computational cost. We introduce a framework for learning when one reduced boundary condition can replace another: paired solutions train a neural network to estimate the resulting domain and boundary errors, and the simpler condition is used only when both predicted errors meet prescribed tolerances. We focus on singular limits in applications, in which a stiff Robin or nonlinear boundary law is replaced by its limiting Dirichlet form. We evaluate the method on a galvanic corrosion problem and other nonlinear stationary and evolution problems.

Abstract: Localized dimensionality reduction improves the scalability of operator learning for high-dimensional partial differential equations (PDEs), but independently decoded local patches can introduce block offsets, interface mismatches, and spurious high-wavenumber content. We introduce Two-Scale Localized PCA-Net, which decomposes the solution into a coarse-global component and local residual corrections. A compact global PCA basis captures domain-scale structure, while nonoverlapping local PCA bases represent the remaining fine-scale residual. A block-balanced latent objective couples the two representations, and optional interface-aware fine-tuning further promotes continuity through reconstruction and trace losses. On Poisson benchmarks, the two-scale representation substantially reduces reconstruction error and visible block artifacts relative to plain and overlap-based localized PCA-Net while approximately halving PCA fitting cost relative to overlap. On heterogeneous Darcy flow, it strongly reduces interface and discrete-residual errors, with more modest reconstruction gains. Ablations show that the primary improvement arises from the two-scale output representation, while interface-aware fine-tuning provides complementary continuity refinement. Overall, separating globally coherent structure from localized residual detail provides an efficient representation for artifact-reduced PDE operator learning.