Papers

This work studies adaptive sparse-grid discontinuous Galerkin (DG) discretizations for the Bhatnagar–Gross–Krook (BGK) model, a kinetic equation posed in four- and six-dimensional phase-space. Standard DG methods are rendered impractical for the BGK model by the curse of dimensionality, motivating compressed representations that adapt to the solution in time. Using the adaptive sparse-grid DG method, we quantify accuracy and compression by comparing the adaptive degrees of freedom to full-grid DG methods and by assessing the resulting kinetic and fluid quantities in both fluid and rarefied regimes. Test cases include a relaxation problem, a multidimensional Sod shock tube, and shear/expansion flows used in prior low-rank BGK studies. To build an efficient Maxwellian evaluation without violating conservation, a central obstacle for structure-perserving BGK simulations, we introduce a hybrid interpolation strategy that exploits velocity separability to recover the correct discrete collision invariants and prove conservation of the resulting discrete collision operator on adaptive sparse grids. Our results show that the adaptive sparse-grid strategy can recover accurate and physically relevant solutions with sharp gradients, and the method reduces the active degrees of freedom by factors ranging from several-fold to several orders of magnitude, with the largest reductions occurring in the six-dimensional examples. All computations are performed with the open-source ASGarD adaptive sparse-grid DG library.

The Solvability Complexity Index (SCI) provides an extensional limit-height formalism for recovering a target map $\Xi$ from finite samples of an evaluation interface $\Lambda$ by finite-height towers of pointwise limits. At first sight this sounds like a typical Type-2 question, so it seems natural to ask how deep and rich the connection of SCI approximation questions and Type-2 computability questions is. This note aims to explain this connection further with explicit examples mostly from the literature.

This paper develops and analyzes an optimal-order fractional backward collocation method for adjoint Volterra integro-differential equations with weakly singular kernels. The backward Volterra structure, together with the weakly singular kernel, induces fractional-power singularities at the terminal endpoint, thereby reducing the classical regularity of the exact solution and causing order deterioration in standard polynomial collocation methods. We first establish a regularity result that characterizes the terminal singular behavior of the solution, showing that the solution is continuously differentiable, whereas its second derivative may exhibit a weak singularity at the terminal endpoint. Motivated by this regularity structure, we introduce a terminally graded mesh and construct a fractional backward collocation scheme whose local approximation space is adapted to the endpoint singularity. Rigorous convergence and superconvergence estimates are derived for both the solution and its derivative. With suitable choices of the fractional parameter and the mesh-grading exponent, the proposed method attains the optimal convergence orders dictated by the local approximation degree. Numerical experiments confirm the theoretical predictions and demonstrate the accuracy and effectiveness of the method for adjoint weakly singular Volterra integro-differential equations.

We develop and analyze a Monte Carlo method for sampling killed anomalous diffusions obtained by time-changing Brownian motion with drift by the inverse of a subordinator. The method targets probabilistic representations of time-nonlocal, including time-fractional, Cauchy–Dirichlet problems on bounded domains. Since inverse subordinators can be sampled exactly in broad classes, while Brownian exit times are generally unavailable in arbitrary domains, we approximate the killed Brownian component by an Euler scheme with discrete boundary detection. We prove a square-root weak error bound with explicit dependence on the Laplace exponent of the subordinator, and derive mean-square and central limit results for the resulting Monte Carlo estimator. A numerical example in the disk and one in a high-dimensional anisotropic shell illustrate the theoretical rates, computation time, and the mesh-free character of the method.

Topology optimization (TopOpt) is a standard tool for structural conceptual design, providing optimal structures with great design freedom. However, it has a high computational costs due to the repeated evaluation of high-fidelity finite element models and frequently produces complex geometries that require extensive interpretation for manufacturability. As an alternative, this work presents a parameterized optimization framework based on Reduced Order Modeling (ROM) and preconditioning techniques. We employ the Empirical Interscale Finite Element Method (EIFEM) coupled with the Discrete Empirical Interpolation Method (DEIM) to construct localized, parameter-dependent reduced operators. From an optimization perspective, this approach can be interpreted as a preconditioning strategy in which inexact gradients are used to accelerate convergence. Furthermore, parameterizing the design space in terms of explicit geometric features, such as inclusion radii or lattice widths, significantly improves the manu facturability of the resulting optimized designs. We assess the performance of the method in terms of computational time, design topology and structural performance with TopOpt as baseline for three differ ent unit cell geometries and three different benchmarks. While the parametric ROM framework requires an initial offline training phase, our comparative analysis demonstrates that it drastically accelerates the online optimization loop and that it can produce even stiffer design in some of our experiments.

We develop an encounter-based formulation of restricted diffusion in a bounded domain with multiple targets, resolving separately the boundary local time accumulated on each target. Unlike the single-target problem, the target projection operators do not generally commute with the governing Dirichlet-to-Neumann (DtN) operator, so that its eigenbasis does not diagonalize the local-time dependence. We overcome this difficulty by decomposing the DtN operator into target-restricted blocks. Multi-dimensional Laplace inversion then yields a convergent switching expansion, whose successive terms describe alternating diffusive transfers between the targets. We also derive an equivalent operator-valued renewal equation. For two targets, diagonalizing the target-restricted blocks provides an explicit spectral representation in terms of two DtN block spectra and an inter-target coupling matrix. When the coupling preserves spectral modes, the switching series can be resummed exactly in terms of modified Bessel functions. In general, as the off-diagonal DtN block is smoothing for separated targets, one can resort to low-rank finite-dimensional approximations. We examine an effective two-mode reduction for small, well-separated targets and identify a regime in which repeated inter-target transfers are progressively suppressed. Probabilistic interpretations and implications for diffusion-controlled reactions are discussed.

We construct variational integrators to study the dynamics of stochastic Hamiltonian systems using their Wong-Zakai approximation. This approach can be used to interpret and justify previous work of Wang et al (L. Wang, J. Hong, R. Scherer, F. Bai. Dynamics and variational Integrators of stochastic Hamiltonian systems. International Journal of Numerical Analysis and Modeling, 6(4), 2009) on stochastic variational integrators. Numerical examples are provided.

Trilinear hexahedral mesh elements are used in finite element analysis to simulate volumetric phenomena. Realism of the simulation mandates that mesh elements through the trilinear map maintain positive volume, or mathematically, that the trilinear map maintains positive Jacobian determinant (jacdet). While jacdet positivity generally needs to be verified in the full element volume, we prove Knupp's conjecture which posits that for trilinear hexahedra, jacdet positivity on the boundary of the element is sufficient to guarantee positivity in the entire element. We further prove that for any valid hexahedron, its globally minimal jacdet must reside on the boundary. Lastly, we prove that the globally minimal jacdet on the boundary of a hexahedron, valid or not, must be achieved within a finite set of candidate points that can be determined by quartic root finding. Combining these results, we can algorithmically validate a hexahedron through evaluation of its jacdet on a finite set of boundary points.

Local characteristic decomposition (LCD) is widely used in high-order numerical methods for hyperbolic systems of conservation laws to reduce spurious oscillations appearing in the computed solutions. The LCD implementation requires a representative average interface state, typically obtained using arithmetic or Roe-type averages of the nearly grid values. We propose an entropy-based LCD (ELCD), in which the states in the two cells adjacent to an interface are considered together with several nearby states in phase space. Among these candidates, we select the state that locally minimizes the entropy and use it as an average interface state for linearizing the flux Jacobian. We incorporate the ELCD into several second-order finite-volume and fifth-order finite-difference schemes. Numerical experiments for the two-dimensional Euler equations of gas dynamics show that the schemes, which utilize the ELCD procedure generally resolve complex wave structures more sharply than their counterparts, which use the arithmetic averages.

The fictitious domain method offers a powerful alternative to circumvent complex mesh generation process by embedding the complex physical domain into a simpler, non-matching background mesh. Among the various strategies to enforce the boundary condition on such a mesh, one possibility is to impose it weakly through a Lagrange multiplier. Existing approaches following this strategy typically require two independent uniform meshes for the domain and the boundary, whose mesh sizes are related by a prescribed ratio. In this work, we relax this requirement by constructing the boundary mesh directly from the trace of the background triangulation. To recover the uniform discrete inf-sup condition under this weaker mesh assumption, we enrich the discrete space of the solution with bubble functions. Our main result is the construction of a suitable operator that allows us to prove this uniform discrete inf-sup condition, thereby establishing the stability of the resulting scheme. As a consequence, we derive optimal a priori error estimates and provide a numerical experiment to validate the theoretical results.

We propose GPU algorithms for high-fidelity simulation of incompressible turbulent flows. Discretization in space is performed with H(div)-conforming high-order Raviart-Thomas finite elements for the velocity and an $L^2$-conforming discontinuous Galerkin approximation for the pressure. In time, a consistent splitting scheme based on higher-order BDF time stepping is used, with convection treated explicitly. In this scheme, a pressure Poisson equation and a symmetric reaction-diffusion-type equation for the velocity need to be solved in each time step. We develop a solution framework with fast matrix-free operator evaluation for all ingredients, combined with multigrid solvers for the Poisson problem, and propose a robust mixed-precision algorithmic framework. A key to mixed-precision efficiency is a least-squares projection to generate accurate initial guesses for the iterative linear solvers, enabling us to work with relative residual tolerances of $10^{-3}$. In this regime, running the solvers entirely in single precision leads to almost no change in overall iteration counts and maintains the crucial turbulence statistics, while showing up to $1.7\times$ speedup over pure double-precision simulations.

An interior penalty nonconforming finite element method is developed for a strain gradient elasticity (SGE) model in arbitrary dimensions, using an $H^1$-nonconforming displacement element built from vector-valued quadratic polynomials enriched by divergence-free face bubbles. Rigorous analysis establishes optimal and robust error estimates with respect to both the Lamé coefficient $\lambda$ and the size parameter $\iota$. The underlying divergence-commuting structure further leads to nonconforming finite element Stokes complexes in two and three dimensions. Numerical experiments support the predicted convergence behavior and parameter robustness, while a circular-hole benchmark under uniaxial tension illustrates the applicability of the method to a multiply connected domain.

We study a finite-domain event for the mirror model on the Manhattan lattice: every boundary-to-boundary trajectory crosses a marked connector an even number of times, while internal cycles are unrestricted. Our main theorem gives a finite-scale confinement criterion: if this event's probability reaches a universal threshold in one even two-square rectangle, then almost surely every trajectory in the plane is periodic. We construct an exactly normalized three-color tensor representation and prove deterministic bounds for approximate contractions using the signed partial-permutation structure of directed continuations. These tools turn numerical evaluations into rigorous finite-volume statements. At mirror density \(p=2/5\), we prove that the boundary-parity probability is nonmonotone in the rectangle width, with a decrease followed by a certified increase. These values remain below the sufficient confinement threshold and do not resolve planar localization at this density. Further experiments track how the measured finite-size minimum varies with density, compare matching-based and entrywise error bounds, and exhibit endpoint-sensitive and cooperative responses to local mirror changes.

Accurate methane source localization and emission-rate estimation at oil and gas facilities require wind fields that resolve site-scale spatial and temporal variability, which cannot be characterized by a single anemometer or coarse operational weather products. We develop a multiscale Weather Research and Forecasting (WRF) framework for a Permian Basin facility that dynamically downscales hourly, 3 km High-Resolution Rapid Refresh (HRRR) fields through four nested domains with grid spacings of 3 km, 1 km, 200 m, and 40 m. The two outer domains use planetary-boundary-layer parameterization, whereas the two inner domains operate in large-eddy-simulation mode. One-minute wind observations from two on-site anemometers are assimilated through wind-only observational nudging in the innermost domain, and the effects of nesting feedback and turbulence-closure choices are examined. Simulations for winter and summer 2025 periods are evaluated using near-surface wind-speed time series, power spectra, spatial fields, and statistical metrics. The baseline simulation reproduces the broad evolution of observed wind events but drifts during weak-wind periods. Observational nudging with one- and two-way nesting reduces the absolute bias of wind-speed by %65 at one sensor and %87 at the other, while reducing the corresponding root-mean-square errors by %33 and %44. Nudging also increases high-frequency energy and resolved spatial gradients, although improvements in turbulence statistics are not uniform. These results demonstrate a practical physics-based pathway for generating high-resolution wind fields for methane-plume modeling while emphasizing the need for independent spatial observations to validate accuracy away from assimilated sensors.

We present an efficient wideband construction of directional $\mathcal{H}^2$-matrices for the electrical field integral equation that, in contrast to existing constructions, supports not only box trees but also geometry-adaptive cluster trees. To accommodate geometry-adaptive cluster trees, we determine the number of directions from the electrical size of each cluster (instead of the level of a box tree), construct the directions using Spherical-Fibonacci points, and establish a hierarchy between the direction sets of a cluster and its children through an angular nearest-neighbor mapping. We construct the nested directional representation using the incomplete adaptive cross approximation, for which we introduce a robustified tree-mimicry pivoting strategy that prevents premature convergence for block-structured matrices arising from certain geometries and meshes. Numerical results demonstrate that the proposed approach achieves the desired accuracy, requires no more storage than the octree-based construction and substantially less when the geometry or discretization is poorly matched to octree clustering, and exhibits the expected $\mathcal{O}(N\log N)$ scaling for high-frequency problems.

The query and key projections $\WQ,\WK$ in attention are almost always trained by Euclidean optimizers with no geometric constraint. We constrain them to the Stiefel manifold and optimize with a Riemannian Adam carrying one scalar second moment per frame—the form of \citet{becigneul2019}, here extended to the compact, non-Hadamard $\St(d,r)$ with a tangent projector, step-norm cap, and polar retraction. Four propositions prove steepest descent in the embedded metric, gradient-scale independence, well-conditioning, and exact $\mathrm{O}(d)$-equivariance. A fifth records that weight decay has identically zero Riemannian gradient on $\St(d,r)$ ($W{=}WI_r$ lies in the normal space), so decay cannot act on the constrained frames. On a CIFAR-10 patch benchmark at $n{=}10\mathrm{k}$ this rule gains $\mathbf{+6.79}$\,pp over AdamW across 12 paired starts ($t{=}38.33$, $12/12$); earlier fixed-step Riemannian SGD gains $+1.97$\,pp, of which $+1.69$\,pp comes from frozen orthonormal initialization alone. The corrected Adam's lead grows with data: $+1.9$\,pp at $n{=}1\mathrm{k}$ to $+6.7$\,pp at $n{=}50\mathrm{k}$. A 12-seed ablation credits all gain to the scale-free step ($+4.63$\,pp, $12/12$), nothing to the projector or equivariance; a targeted $\varepsilon$-sweep causally confirms the mechanism ($-2.6$\,pp at $\varepsilon{=}0.1$, $p{<}0.001$). Two five-seed grokking studies confirm the constrained arm does not grok better than the baseline ($p{=}0.019$, A2 wins): the weight-decay exemption has no grokking consequence. A single-seed pilot exploiting this localization achieves the first stable grokking under slingshot conditions—Stiefel + targeted circuit regularization keeps routing-frame isometry error $10^6\times$ lower than the unconstrained ablation through every collapse.

We propose a deep learning framework for the numerical approximation of weakly singular semilinear elliptic equations subject to homogeneous Dirichlet boundary conditions and positivity constraints. These problems are challenging because the nonlinear source term becomes unbounded as the positive solution approaches zero near the boundary. The proposed approach combines three components. First, the Dirichlet boundary conditions are embedded directly into the neural representation through a hard-constraint formulation. Second, a Softplus transformation preserves positivity in the interior of the computational domain, ensuring that the singular nonlinear term remains well-defined during training. Third, a singularity-aware weighted residual loss emphasizes regions where the predicted solution becomes small and the singular behavior is more pronounced. The method is evaluated in the weak singularity regime through numerical reference solutions and controlled verification tests. Comparisons with the standard residual formulation are used to assess approximation accuracy and convergence behavior. The numerical experiments illustrate the potential of combining hard boundary enforcement, positivity preservation, and singularity-aware residual weighting for weakly singular elliptic problems.

Taming the exponential increase of the Hilbert space dimension with system size in the simulation of gapped quantum lattice models is of the highest relevance for understanding and designing exotic quantum materials, where nonlocal interactions are of particular interest. High-order linked-cluster expansions provide access the solution of the eigenvalue problem for the infinite system, yet rely on the computation of high-dimensional oscillatory lattice sums with a graph structure, only approachable with Monte Carlo methods so far. This work resolves this issue, rendering all required graph lattice sums, referred to as graph zeta functions for kernels involving power-laws, computable. The resulting method reduces the evaluation time for state-of-the art series expansions from tenthousands of core-hours to minutes. After factorizing the lattice sum over blocks, each block is evaluated by the cheapest available strategy depending on its treewidth $\mathrm{tw}$. Basic blocks admit analytic forms in terms of generalized zeta functions. Series-parallel blocks with $\mathrm{tw}\le 2$ can be computed at linear cost in the number of graph nodes and in the size of the momentum grid using a semi-analytical algebra based on Epstein zeta functions and rapidly decaying Fourier series. Finally, for $\mathrm{tw}>2$, the method is combined with tensor-network bucket elimination yielding polynomial scaling of numerical work and memory in momentum grid size with exponents only growing with $\mathrm{tw}$ rather than with the number of vertices. Through use of FFT, the full momentum grid is recovered at the cost of a single momentum evaluation. We provide a detailed analysis of the precision and runtime of our method against analytic and numerical benchmarks. We further reproduce published Monte Carlo data for the transverse-field Ising model on different 1D, 2D, and 3D lattices, obtaining full agreement.

We propose a two-dimensional benchmark for fluid-structure-contact interaction consisting of a deformable elastic disk falling under gravity within a viscous incompressible fluid and rebounding in the vicinity of the bottom wall. Solid deformability is essential, as rigid solids do not rebound in this framework. Besides this, the setting is deliberately kept simple to facilitate reproduction. The configuration is particularly challenging due to the well-known no-contact paradox, which can lead to a contactless rebound and forces numerical methods to resolve a vanishingly thin fluid layer in the near-contact region, making the dynamics highly sensitive to the spatial and temporal discretizations. In addition to no-slip boundary and interface conditions, a reduced porous modeling of surface roughness, either on the disk boundary or on the bottom wall, is also considered; this circumvents the no-contact paradox and enables genuine contact. An energy balance law is derived theoretically for all three cases. Eight numerical methodologies, developed by five research groups and spanning different model formulations, numerical methods, and codes (including both fitted and unfitted discretizations), are applied to the benchmark at several levels of spatial and temporal refinement. Quantities of interest of varying complexity are collected and compared, showing close agreement during the falling phase and increased sensitivity in the near-contact and rebound regimes. The setting and the results provide a suitable reference for the systematic assessment of fluid-structure-contact interaction solvers. The time histories of all quantities of interest for every approach and refinement level are provided as supplementary material.

In this paper, we introduce a novel multi-fluid approximation scheme for the one-dimensional Vlasov-Poisson system, that is relevant for arbitrarily long time intervals. For low regularity solutions, we establish an error estimate in the Wasserstein-1 distance, obtaining a convergence rate of order $O(h^{\frac{2}{3}})$, where h denotes the velocity grid size. This establishes the consistency, as the grid size tends to zero, between the pressureless Euler-Poisson and the Vlasov-Poisson systems. Numerical illustrations are given to illustrate the efficiency of our approach.

We present a fixed nondegenerate tetrahedron with a nested largest-dihedral-anglebisection branch that collapses onto a unit edge. The selected dihedral angles to bisect are uniquely largest ones at every bisection step. Some face angles of the branch elements tend to zero, the corresponding diameter-to-inradius ratio diverges, and the diameters of the elements produced do not tend to zero, although every dihedral angle stays uniformly away from both 0 and {\pi} and all face and dihedral angles satisfy a uniform maximum-angle bound. The proof is based on an exact invariant-range argument. This example shows that the largest-dihedral-angle bisection algorithm can produce degenerating tetrahedral partitions.

In this letter, a novel 3-D Finite-Difference Time-Domain (FDTD) framework is proposed that circumvents the costly volumetric scaling of conventional FDTD methods. By representing the electromagnetic fields as low-rank Quantized Tensor Trains (QTT), the memory requirements scale logarithmically with system size. Moreover, the various numerical operations that constitute the FDTD scheme can be efficiently implemented in this format, with their computational cost also exhibiting a logarithmic complexity. The simulation of systems in open space is enabled by the inclusion of a uniaxial PML. A validation example demonstrates that the method achieves excellent accuracy compared to the traditional full-grid (FG) FDTD method, while significantly reducing the required computational resources. Memory savings of several orders of magnitude are obtained, highlighting the potential of the proposed framework for the simulation of large multiscale electromagnetic systems.

Common finite element discretizations for the incompressible Stokes equations with a continuous pressure approximation – the MINI element, the Taylor–Hood element, or stabilized equal-order elements – are not pressure-robust: the velocity error is polluted by the pressure best-approximation error, multiplied by the inverse viscosity. The established remedy replaces the test function in the momentum equation by a divergence-preserving reconstruction operator; for continuous pressure approximations this construction is cumbersome, requiring the solution of local problems (Lederer, Linke, Merdon & Schöberl, SIAM J. Numer. Anal., 2017). We propose a simpler alternative: interpolating the force by a Nédélec interpolation operator of suitable order instead of reconstructing the test function. Exploiting the commuting diagram property linking Nédélec and Lagrange interpolation operators, we show that the resulting method is pressure-robust for any $H^1$-conforming, continuous-pressure discretization. We give the complete analysis for the MINI element, and then for stabilized equal-order $P^kP^k$ elements of arbitrary degree $k$, proving that the interpolation-induced consistency error is pressure-robust and converges of optimal order. We validate the theory by several numerical examples.

We develop and analyze a fully discrete, globally divergence-free hybridizable discontinuous Galerkin (HDG) method for a gradient-based Smagorinsky model. The method combines backward Euler time stepping, interior-penalty discretizations of molecular and nonlinear eddy diffusion, and an upwind convective flux. The discrete velocity is $H(\operatorname{div})$-conforming and pointwise divergence-free, which yields pressure robustness. For sufficiently large penalty parameters, we prove unconditional energy stability and existence of a discrete solution, and establish uniqueness under additional smallness conditions. A velocity error estimate is derived without explicit inverse powers of the molecular viscosity. The nonlinear facet residuals are controlled using local trace-approximation estimates and a viscosity-independent facet penalty. We retain the dependence of the discrete Gronwall factor on the filter scale and the mesh size; a mesh-uniform bound follows under suitable solution regularity, a fixed time-step margin, and the scaling $\delta=O(h)$ on quasi-uniform meshes. The reported manufactured-solution results are consistent with the resulting pre-asymptotic error bounds. Further flow examples illustrate the dissipative behavior of the method and are distinguished from the boundary conditions and parameter range covered by the analysis.

This paper aims to develop efficient numerical approximations for a class of stochastic differential equations with state-dependent fast switching processes. The direct Euler–Maruyama (EM) scheme fails when the scaling parameter is small. Based on the heterogeneous multiscale method of \cite{e2005analysis}, we propose three algorithms and prove their strong $L^p$-convergence with explicit rates for any $p\geq 2$. In the first algorithm, we combine the averaging principle with an EM scheme for the averaged equation, where the invariant measure of the Markov chain can be explicitly obtained by solving a linear system. However, computing this invariant measure incurs cubic cost as the number of switching states increases. To avoid solving large linear systems, we approximate the invariant measure instead. Thus the second and third algorithms both combine a macroscopic EM scheme for a modified averaged equation with micro-solvers that estimate the averaged drift. More precisely, in the second algorithm, a discrete-time Markov chain is simulated with a micro time step, and the averaged drift is obtained by averaging over finitely many micro transitions. However, both the first and second algorithms only work when the switching process has finite states. Therefore, we introduce a third algorithm that allows for switching processes with countably infinite states, in which exact continuous-time Markov chain is generated via the Gillespie algorithm, and the averaged drift is computed by exact time averaging over a specified interval. Numerical experiments verify the theoretical results and demonstrate the computational advantages of these three algorithms.

The classical Babuška paradox shows that solutions of simply supported plate problems on polygonal approximations of a curved domain may converge to an unintended limit. We develop a boundary-corrected $H({\rm divdiv};\mathbb S)$-$L^2$ mixed finite element method for the simply supported Kirchhoff-Love plate problem. The correction uses only the discrete boundary geometry. Introducing the bending moment as an independent unknown allows the condition $M_{nn}=0$ to be imposed directly, while an edgewise mean constraint on the effective shear suppresses the leading geometric inconsistency. This constraint improves the boundary consistency error from $\mathcal O(h^{1/2})$ to $\mathcal O(h^{3/2})$. The associated boundary corrections act in the kernel of $\rm{divdiv}$, leaving the discrete equilibrium equation unchanged and yielding a uniformly stable scheme. Under suitable regularity assumptions, we prove $L^2$-error estimates of order $h^{3/2}$ for the bending moment and the broken Hessian of the postprocessed displacement, and of order $h^2$ for the displacement. The analysis covers multiply connected domains and polygonal approximations whose boundaries may cross the physical boundary. Numerical experiments confirm these rates and the improvement over the uncorrected method.

Latent neural operators improve the efficiency of operator learning for partial differential equations (PDEs) by performing the main computation on compact latent representations. However, directly compressing the input representation to obtain such compact representations may discard solution-relevant spatial information, especially for PDE solutions with multiscale structures. To address this problem, we propose HiLNO, a hierarchical latent neural operator that constructs a fine-to-coarse-to-fine latent space and further introduces multi-scale supervision (MSS) and anisotropic Gaussian attention. The hierarchy mitigates potential information loss during compression, while MSS aligns intermediate predictions with downsampled target fields, encouraging solution-relevant structures to be captured across multiple spatial scales. Anisotropic Gaussian attention enables feature transfer across the hierarchy, making HiLNO applicable to general geometries. Experiments on representative PDE benchmarks and a large-scale automotive aerodynamics task show that HiLNO achieves competitive predictive accuracy, while reducing the parameter count by an average of 84.4% and FLOPs by an average of 69.2% compared with LinearNO. Additional experiments demonstrate effective generalization to unseen spatial resolutions. Code is available at https://github.com/JcLimath/HiLNO.

Phase-field trajectories governed by the Allen–Cahn equation determine the effective force $G=F'/\eps^2$, but not the free-energy potential $F$ and interfacial scale $\eps$ separately. We introduce surface-tension-calibrated Allen–Cahn identification (STAC), a two-step method that resolves this ambiguity using an independently supplied surface tension. First, space–time weak moments produce a linear inverse problem without differentiating the measured field. A Bernstein representation with positive coefficients enforces two stable pure phases and a single central energy maximum. Second, the surface tension fixes the remaining scale algebraically. We prove scale equivalence of the forward model, uniqueness of the calibrated pair within this class, and finite perturbation bounds that remain valid when the potential vanishes at the pure phases. Experiments with resolved synthetic phase fields separate the effects of the structural constraint, ridge regularization, and constitutive approximation error. Across 80 shared noisy records, unconstrained degree-five regression produces six potentials with an incorrect well structure. Positivity combined with ridge regularization eliminates these violations and reduces the mean force error from 8.70% to 5.66%. Lower-degree models are nevertheless more accurate in some tests, and an admissible law outside the positive Bernstein family shows how the constraint can introduce approximation bias. STAC thus recovers a structurally admissible, physically scaled potential while making explicit which information comes from the observed dynamics and which comes from the independent surface-tension datum.

We introduce a new finite element discretization for a mixed formulation of the two-dimensional incompressible Stokes equations with symmetric viscous stresses. The method is based on a mass-conserving stress-yielding formulation (Gopalakrishnan, Lederer, Schöberl; A mass conserving mixed stress formulation for the Stokes equations; IMA J. Numer. Anal., Vol 40, 2020) with symmetric stresses, with the key novelty that the symmetry of the stress tensor is enforced exactly as an intrinsic property of the discrete space. This space is constructed on the Clough–Tocher macroelement split and consists of matrix-valued functions whose normal-tangential components across macroelement interfaces are continuous. Since the exact symmetry is enforced via local polynomial bubble functions on the subelements, which can be eliminated via static condensation, the number of globally coupled degrees of freedom coincide with those of the original mass-conserving stress-yielding formulation. The resulting pressure-robust method is stable without requiring any additional enrichment. We establish optimal convergence rates for all the variables and present numerical experiments that confirm the theoretical results.

It is mathematically proved that, for a grid under a boundary-preserving piecewise-smooth deformation, the positivity of the Jacobian determinant in each closed element is also a sufficient condition for the deformed grid to remain valid. Most of the proof is devoted to showing that the element-wise positivity of the Jacobian indeed induces a local isomorphism in the grid-deformation setting. This result provides users of ALE methods and of the generalized Lagrangian method with a clean answer: checking the Jacobian alone is enough to determine the validity of the deformed grid.