Papers

We introduce a Hybrid High-Order framework for differential forms on general polyhedral meshes and characterise stability-admissible face spaces in terms of the traces of the polynomial kernel of the exterior derivative. The analysis is based on a cellular-to-hybrid transfer mechanism in which the global part of the stability problem is reduced to a single cellular cochain Poincaré problem, and the resulting control is propagated to the hybrid level using stable polynomial skeletons and local polynomial completions. This mechanism yields a uniform stable conforming lifting and a uniform discrete Poincaré inequality for every form degree in arbitrary space dimension. The lifting preserves the prescribed reconstructed exterior derivative after projection, whereas the Poincaré inequality controls the broken cell polynomial form modulo the continuous conforming kernel of the exterior derivative. Both results hold on domains with arbitrary topology and uniformly over stability-admissible choices of the face spaces. In three dimensions, the conforming-kernel estimate recovers the hybrid Poincaré–Wirtinger inequality for the gradient, the second hybrid Weber estimate under the standard gauge, and the corresponding divergence-kernel estimate.

Fractional diffusion models provide a powerful framework for describing anomalous transport phenomena in heterogeneous porous media. The Mobile-Immobile model is a fundamental approach for characterizing such anomalous diffusion, specifically addressing the delayed solute transport caused by mass transfer between mobile and immobile regions. In this study, we develop and analyze two fully discrete numerical schemes for the one-dimensional time-fractional Mobile-Immobile model governed by the Caputo derivative of order $\alpha\in (0,1)$. The spatial discretization is carried out using a non-symmetric interior penalty discontinuous Galerkin method, while the temporal derivative is approximated by the Crank-Nicolson L1 and L2-$1_\sigma$ formulas, resulting in two distinct schemes with high-order accuracy. Rigorous stability and error analyses are established, showing that the methods achieve optimal convergence rates depending on the regularity of the exact solution. Numerical experiments verify the theoretical predictions and demonstrate the efficiency and robustness of the proposed schemes. The presented framework provides a reliable and accurate approach for simulating complex fractional transport processes in porous media and related fields.

In this paper, we propose an efficient and robust fifth-order finite-volume Hermite weighted essentially non-oscillatory (HWENO) scheme with gradient reconstruction for the compressible Navier–Stokes equations. The key idea is to construct weak-derivative moments for each dissipative variable by taking the arithmetic average of two HWENO interface traces, which are then processed componentwise through the standard scalar nonlinear HWENO reconstruction formula, with only the input moments differing from those used in the convective reconstruction, thereby supplying the gradients needed for the viscous fluxes. This reconstruction achieves fifth-order accuracy using only candidate polynomials of degree at most four, avoids the order reduction that would result from direct differentiation, and introduces negligible additional computational cost. To ensure robustness, especially in extreme test cases involving strong shocks or low densities, we apply a positivity-preserving limiter to the conservative states while retaining the nonlinear HWENO reconstruction for the viscous gradients. This combination is essential for stable computations in flows with strong gradients. Both the conservative and gradient reconstructions share the same compact stencils and candidate-polynomial structure, enabling a single implementation routine without introducing additional algorithmic complexity. Numerical results confirm that the proposed scheme delivers fifth-order accuracy and high resolution, offers competitive computational efficiency, and robustly resolves challenging finite-Reynolds-number Navier–Stokes flows in which the positivity-preserving limiter is actively engaged.

This article proposes and analyses a quadratic $C^0$ interior penalty method for the displacement obstacle problem of the von Kármán plate. The discrete space consists of Lagrange $P_2$ finite elements and the obstacle constraint is imposed at the vertices. The trilinear form of the von Kármán bracket is modified by terms on the edges so that it is bounded in the discrete energy norm. The well-posedness of the discrete problem, namely the existence of a discrete solution and its uniqueness under a smallness condition on the data, is established. The Sobolev and Friedrichs constants of the discrete energy norm are quantified with an explicit dependence on the mesh size. The main result is an error estimate of order $\mathcal{O}(h^{\alpha})$ in the discrete energy norm, where $1/2<\alpha\le1$ is the index of elliptic regularity of the biharmonic operator on the polygonal domain. Numerical experiments on a square and on an L-shaped domain confirm the predicted rates. The coincidence set has positive measure in one example and empty interior in another. The experiments also show how the penalty parameter affects the rates and identify a threshold in the size of the obstacle beyond which the iterative solver fails on fine meshes.

In this paper, we develop a high-order, energy-stable and charge-conservative Lagrangian finite element method for variable-density incompressible inductionless magnetohydrodynamic (MHD) equations. The method utilizes the moving high-order curved tetrahedral mesh to track the material interface. Second-order Backward Differentiation Formula (BDF2) is used for the temporal discretization of material derivative, together with the second-order Adams–Bashforth method (AB2) for the update of control points of the meshes. High-order isoparametric Taylor-Hood elements with grad-div stabilization are used for the velocity-pressure pair. While, to ensure the discrete charge conservation, high-order parametric $\BH(\Div)$-conforming element is adopted for the current density. In the absence of external force, the unconditional energy-stability of the fully discrete scheme is proven. Finally, 3D numerical experiments are conducted to confirm the expected high-order accuracy for smooth solutions, the energy stability property and the capability of the proposed method.

In this paper, we present a fifth-order finite difference divergence-free Hermite weighted essentially non-oscillatory (HWENO) scheme for the ideal magnetohydrodynamics (MHD) equations. In this framework, both the solution and its spatial partial derivatives are evolved in time and jointly employed in the spatial reconstruction procedure. A major challenge in MHD simulations is preserving the divergence-free constraint of the magnetic field, which is generally violated by standard numerical methods designed solely for hyperbolic conservation laws. To address this issue, we first solve the MHD equations within the HWENO framework for hyperbolic conservation laws, yielding a magnetic field divergence that remains zero up to high-order accuracy in smooth regions. Subsequently, we apply a correction that evenly distributes the divergence error among the partial derivatives involved in the divergence-free constraint, thereby rendering the magnetic field discretely divergence-free at the new time level. This approach offers several advantages. First, the scheme retains the conservation property, as only the partial derivatives of the numerical solution are corrected, leaving the conserved variables unchanged. Second, the correction applied to the partial derivatives of the magnetic field components introduces only a high-order perturbation, thereby preserving the overall accuracy. Third, the divergence-free treatment significantly enhances robustness, as most benchmark test cases cannot be run stably without such a correction. Fourth, the correction is a simple linear operation applied at each stage of the time integration, incurring negligible additional computational cost. Extensive numerical experiments demonstrate the accuracy, resolution, efficiency, effectiveness, and robustness of the proposed scheme.

In recent years, neural networks have significantly advanced numerical solutions of partial differential equations (PDEs). However, solving PDEs with discontinuous solutions, such as hyperbolic conservation laws, remains challenging for neural network-based methods such as physics-informed neural networks (PINNs). Existing methods often rely on strong prior assumptions such as knowledge of discontinuity locations, or they introduce artificial smoothing terms that degrade accuracy. However, accurately solving these conservation laws and predicting the formation and propagation of discontinuities in solutions is crucial in many practical applications, including gas dynamics and traffic flow modeling. In this paper, we introduce a novel Weak-Entropy PINN (WEPINN) framework for hyperbolic conservation laws with discontinuous solutions. The method enforces the governing equations in their weak (integral) formulation and incorporates the entropy condition to select the physically admissible solution, while employing the discrete fast Fourier transform (DFFT) for efficient numerical integration. Our method is tested through extensive numerical experiments on a variety of scalar conservation laws and systems of conservation laws in one and two dimensional spaces. These experiments demonstrate that our method can accurately resolve sharp discontinuities while effectively capturing interactions between multiple shock and rarefaction waves.

A thin shell carries load through its shape and, at finite deflections, its stiffness changes with the load itself, so a shape optimum found with a linear model can be far from optimal. We present an automated framework that embeds the fully geometrically nonlinear shell response in the shape-optimisation loop. The forward model, a five-parameter nonlinear Naghdi shell stabilised against locking by partial selective reduced integration, operates directly on a discrete (faceted) triangulation with a numerically recovered director field – dispensing with the exact mid-surface parameterisation of isogeometric approaches, a chart that ceases to exist once the geometry itself is the design variable. Implemented in Firedrake, the model generates its residual, consistent tangent and adjoint automatically; shape derivatives, computed by algorithmic differentiation through the full load-continuation solve, drive the Fireshape/ROL trust-region optimiser. The forward solver reproduces the Sze/Abaqus benchmark for a clamped semi-cylindrical shell under a point load, capturing the progressive stiffening that a geometrically linear model cannot reproduce. The optimisation is validated against the COMSOL benchmark, a sheet-metal bracket under bending: the framework develops the same off-mid-plane corrugation mechanism and attains an 87% reduction of elastic strain energy within the prescribed displacement budget, matching the benchmark's magnitude and area change. Applied to the curved semi-cylinder, it forms a smooth stiffening crease that reduces the shell's average deflection under load by 78%.

High-order shock-capturing schemes often face severe computational bottlenecks due to wide stencils, stagewise nonlinear weights, and expensive local characteristic decompositions. To address this, we propose a compact finite-average-moment (FAM) framework for hyperbolic conservation laws that completely decouples the formal spatial order from the number of evolved local degrees of freedom. By evolving only a $\mathbb{P}^1$ moment state (i.e., cell averages and scaled first-order moments), we achieve up to sixth-order accuracy ($k=3, 4, 5, 6$) via compact, polynomially exact linear moment reconstructions. A key algorithmic innovation is consolidating the nonlinear stabilization into a single, derivative-free oscillation-elimination (OE) procedure applied only after the final Runge–Kutta stage. This OE procedure exactly preserves cell averages while applying an explicit exponential correction to the first-order moments, avoiding repeated stagewise nonlinear weights or limiting. Fourier analysis of the linear backbone demonstrates $k$th-order accuracy, strong spurious mode damping, and $(k+1)$th-order cell-average superconvergence. Extensive numerical experiments on scalar laws and the Euler equations verify the expected smooth accuracy and robust, nonoscillatory shock resolution. Notably, by performing componentwise reconstruction directly in conservative variables and streamlining stabilization, the FAM schemes deliver significantly lower complete-run wall-clock times compared to the multi-resolution weighted essentially non-oscillatory (MR-WENO), unified-stencil Hermite WENO (HWENO-U), and OE-HWENO methods.

The finite element immersed boundary method (FE IBM) with a distributed Lagrange multiplier is an attractive framework for fluid structure interaction (FSI) because it couples an Eulerian description of an incompressible fluid to a Lagrangian description of an immersed solid on two independent, nonconforming meshes, avoiding the costly remeshing required by body fitted arbitrary Lagrangian Eulerian methods. When the immersed solid is modelled as a fully incompressible hyperelastic material, however, a direct finite element discretization of the deviatoric solid stress fails to enforce the Lagrangian incompressibility constraint, and the computed structure exhibits spurious volumetric instabilities and locking. In this work we present a mixed formulation of the distributed Lagrange multiplier FEIBM that restores volumetric stability. Following the theory of nearly incompressible hyperelasticity, we augment the solid stress with a volumetric contribution derived from a dilatational strain energy and introduce an additional solid pressure field, enforced weakly, that plays the role of the Lagrange multiplier for the Lagrangian incompressibility constraint J = 1. The resulting formulation is discretized in space by finite elements and in time by an unconditionally stable semi implicit scheme, and is implemented as a reusable Immersed Boundary physics kernel within the GRINS multiphysics framework, built on the libMesh finite element library. The method is verified on three FSI benchmarks, an elliptically displaced thick ring returning to equilibrium, a radially stretched incompressible ring, and a disk falling under gravity in a viscous fluid, for which the mixed formulation removes the volumetric failure observed with the unstabilized formulation and reproduces the analytical terminal velocity of the falling disk to within 1%.

Ensuring the bijectivity of spline-based parameterizations is fundamental in geometric modeling and isogeometric analysis, as invalid mappings may lead to self-intersections, singular Jacobians, and numerical instability. While T-splines offer enhanced flexibility through local refinement, this flexibility also makes bijectivity verification significantly more challenging. In this work, we propose a rigorous and efficient framework for bijectivity analysis of rational T-spline surfaces based on Bézier extraction. The key idea is to reformulate the T-spline representation into a collection of element-wise rational Bézier patches, on which the Gram determinant of the mapping admits a Bernstein polynomial representation. This enables a coefficient-based analysis of local regularity by exploiting the convex hull and positivity properties of the Bernstein basis. Based on this formulation, we derive a sufficient condition for bijectivity from the nonnegativity of Bernstein coefficients, together with a necessary condition based on the sign consistency of corner coefficients. For cases where these conditions are inconclusive, we introduce a hierarchical subdivision strategy that progressively localizes ambiguous regions and resolves them through refinement. The proposed method provides a certified and adaptive procedure for bijectivity verification that avoids dense numerical sampling and remains computationally efficient. Numerical experiments on complex T-spline geometries demonstrate that the approach accurately detects both valid and near-degenerate configurations, while scaling effectively to large models with thousands of rational Bézier patches. The framework is fully compatible with standard isogeometric analysis workflows.

We analyze a Galerkin spectral method applied to the Vlasov-Poisson (VP) system based on time-independent asymmetrically-weighted (AW) Hermite functions in velocity. The VP system is written as an hyperbolic system using AW Hermite functions in velocity. Unlike the classical Petrov-Galerkin spectral method, which is not stable, the proposed Galerkin spectral method admits a natural stability in the unweighted L2 norm. This stability property enables a rigorous convergence analysis of the method. For sufficiently regular solutions with exponential decay in velocity, we establish error estimates between the exact and numerical solutions and prove convergence of the Galerkin spectral method. In particular, the method achieves spectral convergence in Sobolev spaces, with convergence rates determined by the regularity of the exact solution. Since the Gram matrix arising from the proposed Galerkin method is dense, we derive an equivalent form that retains the same approximation properties while preserving the sparsity structure similar to the classical Petrov-Galerkin method.

Computing (adjoint) sensitivities of flows with shocks is a longstanding problem in computational fluid dynamics. The spurious sensitivities due to shock sensors and limiters frequently force practitioners to accept the errors incurred by "freezing" limiters and shock sensors. The recently proposed information geometric regularization (IGR) is an inviscid, PDE-based regularization of the compressible Euler equations that replaces shocks with smooth profiles without damping fine-scale structures. This work derives the forward and adjoint sensitivities for the IGR system with periodic boundary conditions and demonstrates their convergence, under grid refinement, to the sensitivities obtained by finite differences or automatic differentiation through the forward solve.

The stabilising effect of multiplicative noise for stochastic differential equations, though counterintuitive, has been observed and investigated extensively in last decades. In practice, it is desirable to know if such stabilisation holds also for the discretised setting. In this paper, we address this problem by means of sharp upper and lower estimates for Lyapunov exponents of Milstein schemes for non-autonomous stochastic differential equations. These estimates provide precise large time behaviour in both almost sure and $p$-moment sense. In particular, our results show the preservation of stabilisation from the continuum setting to the discretised setting. One main idea of our analysis is to exploit the second order term concerning the stochastic noise from the Milstein scheme to obtain precise estimates for Taylor expansions of logarithmic and power functions.

Recent research has demonstrated the possibility of exploiting the thermodynamics of coupled electrical oscillators to implement computational tasks such as matrix inversion. While physical implementations rely on thermal noise to drive equilibration, we show that the underlying dynamics reduce to a deterministic iterative algorithm. Building on the framework of Aifer et al., we analyze the moment evolution of the Ornstein-Uhlenbeck process governing thermodynamic symmetric positive definite (SPD) matrix inversion. We prove that to a first-order approximation, the covariance dynamics are mathematically identical to preconditioned gradient descent on the Frobenius norm of the residual $\tilde{A}^{-1}A-I$. This equivalence demonstrates that thermal fluctuations, while essential for physical thermodynamic hardware, are algorithmically redundant for convex problems with a single global minimum. We validate the resulting algorithm against Thermox (Duffield et al.), a stochastic thermodynamic simulator, achieving speedups exceeding 100,000-fold while remaining competitive with the Newton-Schulz iteration. We also demonstrate acceleration through Schur complement techniques. These results establish a rigorous link between analog thermodynamic computing, statistical physics, and deterministic optimization methods.

We study an inverse initial-data problem for a quasilinear transport equation with nonlinear and memory effects. The unknown initial state is reconstructed from time-dependent measurements on the outflow boundary, with prescribed inflow data. We first apply a Legendre–exponential time-dimensional reduction to transform the governing equation into a finite nonlinear system in space. We then develop a Carleman-weighted and Tikhonov-regularized Picard method. At each iteration, the nonlinear terms are evaluated using the previous iterate, leading to a linear minimization problem with a unique solution. A Carleman estimate for the principal transport operator is used to prove that, for a sufficiently large Carleman parameter, the resulting Picard map is contractive on a prescribed admissible set. Consequently, the method converges from an arbitrary initial guess in that set. We also establish stability with respect to noisy outflow data. These analytical results concern the truncated and regularized reduced problem. Two-dimensional numerical experiments demonstrate accurate reconstruction of single and multiple inclusions, robustness with respect to noise, and rapid convergence of the Picard iteration.

In this paper, we study kernel ridge regression using the data-driven diffusion maps (DM) kernel, which is constructed through algebraic manipulations of the diffusion maps algorithm. Under the assumptions that the data lie on a manifold and are sampled uniformly, we prove that the appropriately scaled DM kernel converges uniformly to the heat kernel on the manifold for sufficiently large times as the dataset size increases and the kernel bandwidth is scaled appropriately. Consequently, the limiting Reproducing Kernel Hilbert Space (RKHS) induced by the DM kernel coincides with the RKHS associated with the heat kernel on the manifold. We further show that the RKHS induced by the DM kernel is isometrically isomorphic to the RKHS of the Gaussian kernel, which is continuously embedded in the RKHS of a Matérn kernel whose norm is equivalent to an appropriate Sobolev norm. This result implies that standard risk bounds for kernel ridge regression applicable to Matérn kernels also apply to the DM kernel. Finally, we provide numerical results that (1) validate the convergence of the heat kernel approximation, (2) demonstrate the greater expressiveness of the DM kernel compared to the Gaussian kernel for supervised learning over a larger class of functions on manifolds with boundary, and (3) demonstrate the advantage of the DM kernel over the Gaussian kernel in learning functions with varying frequencies and co-dimensions.

The discretization of second-order (viscous) terms in Discontinuous Spectral Element Methods (DSEMs) typically relies on an auxiliary gradient variable, whose treatment at element interfaces affects the accuracy and stability of the scheme. The Bassi-Rebay (BR1) formulation is attractive for its simplicity and parameter-free character, but suffers from sub-optimal convergence at even polynomial orders and requires an extended five-element stencil. Inspired by Huynh's Flux Reconstruction formulation, we develop a compact, fully-centered scheme for second-order fluxes within the Spectral Difference (SD) method. The proposed approach modifies the reconstruction of the auxiliary gradient using interface-dependent, one-sided continuous fluxes, reducing the stencil from five to three elements while preserving the centered and parameter-free nature of BR1. The formulation is developed in one dimension and extended to multiple dimensions. Temporal eigenanalysis is used to characterize its dissipation and dispersion properties, including the effects of interior penalty terms. Numerical tests consider the linear diffusion equation, an under-resolved localized Dirac's delta, the nonlinear porous medium equation, and implicit large-eddy simulations of the three-dimensional Taylor-Green vortex at $\mathrm{Re}=1600$ and $5000$. The compact scheme restores the expected convergence order for all polynomial degrees, including even orders, and reduces spurious oscillations in under-resolved and nonlinear regimes. It also remains stable in turbulent cases where the standard formulation fails, owing to improved damping of high-wavenumber numerical modes. The proposed approach provides an attractive compact alternative to BR1 for second-order fluxes in the SD method.

This paper introduces a novel block preconditioner for a class of double saddle point systems arising from the finite element discretization of liquid crystal director modeling. We first establish sufficient conditions to guarantee the convergence of the underlying iterative scheme. A spectral analysis of the preconditioned matrix is then performed, leading to an effective strategy for parameter selection. The new preconditioner is used to accelerate the convergence of the flexible version of GMRES. Experimental results confirm the strong performance of the method.

We consider smooth convex minimization over the spectrahedron using Frank-Wolfe-type methods based only on extreme-eigenvector computations. In our recent work \cite{garber2026randomized} we presented the first ambient-dimension-independent linear convergence rate under quadratic growth. However, the method makes an additional strong strict complementarity assumption, it is randomized, its linear rate holds only after a burn-in phase and in expectation, and it requires the objective smoothness constant. We show that these limitations can be removed. Assuming quadratic growth and that all optimal solutions have the same rank, but without assuming strict complementarity, we give a deterministic and parameter-free Frank-Wolfe-type method with a global ambient-dimension-independent linear convergence rate.

How expressive is prompting a transformer? Answering this question is important for separating the roles of prompting, architecture, and pretraining in transformer models, and for determining whether task-specific behavior must be stored in model weights or can instead be induced at inference time through the prompt. We show, in an approximation-theoretic sense, that pretraining is optional: a single-layer softmax attention network with random, untrained weights can approximate any Hölder function on a compact manifold when steered by an appropriate soft prompt. Guided by the connection between softmax attention and kernel methods, we construct explicit soft prompts (a prompt per target function, independent of the query) as solutions to linear systems matching attention logits to Gaussian kernel exponents, under which the frozen transformer emulates the classical Nadaraya-Watson kernel estimator. The construction requires only a mild rank condition on the weights, which we show holds almost surely under Gaussian initialization. The prompted network inherits the theoretical guarantees of kernel regression, leading to universal approximation theorems with minimax-optimal rates that depend on the intrinsic dimension. We further quantify the cost of prompting, exposing a tradeoff between the norm of the constructed soft prompt tokens, prompt length, and hidden dimension. Numerical experiments corroborate the constructions and predicted rates.

We revisit the nonlinear two-scale Dynamic Diffusion (DD) finite element formulation mathematically analyzed by Santos et al. (2021) for stationary advection–diffusion–reaction problems. We establish an explicit local Lipschitz estimate for the artificial diffusivity, with a constant of order $h_T$, and use it to prove uniqueness of the discrete solution for sufficiently fine meshes. By separating the approximation error in the energy norm from the contribution associated with artificial diffusion, we derive an optimal first-order a priori energy-norm estimate for continuous piecewise linear finite elements enriched with simplex bubble functions. We further prove that the square root of the nonlinear artificial dissipation is $O(h)$. Consequently, the combined energy-error and artificial-dissipation measure also converges with first order, sharpening the previously available $O(h^{1/2})$ estimate. A smooth manufactured problem, considered with and without reaction, corroborates the predicted convergence rates. A second, strongly advection-dominated problem with sharp outflow layers illustrates the stabilizing behavior of the DD formulation.

In this paper we provide Monte Carlo and deep neural network approximations for stochastic representations of solutions to linear elliptic partial differential equations with constant diffusion, drift and killing. Building on the modified Walk-on-Spheres algorithm of Beznea et al. (arXiv:2209.01432), we introduce Monte Carlo estimators that explicitly incorporate sampled random times arising in the analyzed stochastic representations. We establish uniform error bounds for these estimators and show that, under suitable assumptions, a prescribed approximation accuracy is achieved with sample complexities growing at most polynomially in both the inverse accuracy and the problem dimension. Furthermore, we prove a deep neural network approximation result for the stochastic representations. Assuming suitable neural network representations of the boundary data and the distance function to the boundary, we use the constructed Monte Carlo to design deep neural networks that approximate the representation uniformly with a number of parameters growing at most polynomially in the inverse accuracy and the problem dimension. These results extend previous complexity analyses to a broader class of elliptic equations involving drift and killing.

Recently, tensor-based physics-informed neural networks (T-PINNs) have received increasing attention. However, existing T-PINNs still face a fundamental challenge: they mainly rely on pre-specified low-rank tensor decompositions with manually tuned ranks, which limits their ability to capture the underlying structures of multivariate solution functions and hinders their practical deployment. To address this challenge, we propose a hierarchical rank-evolving (abbreviated as HRE) representation for multivariate functions, which endows us to faithfully capture the underlying structure of the targeted multivariate function accompanying with automatic rank determination. Concretely, in the hierarchical design of HRE representation, the target multivariate function is decomposed as a small-scale inner tensor with a set of univariate functions along each mode, where a customized tensor network decomposition can be readily deployed to capture the underlying structure of the small-scale inner tensor. In HRE representation, the crucial hyperparameters, ranks, can be adaptively revealed during the decomposition, freeing us from manual rank tuning and making HRE practically applicable to real-world problems. Besides, we build the HRE-PINNs correspondingly. Extensive numerical experiments, including high-dimensional static problems (Helmholtz equation and Poisson equation), nonlinear time-dependent problems (Klein-Gordon equation), and complex fluid-dynamics problems (flow mixing equation and Navier-Stokes equation), demonstrate that HRE-PINNs consistently outperform existing state-of-the-art approaches in terms of accuracy.

We develop an unfitted ghost finite element method for elliptic interface problems with discontinuous diffusion coefficients and apply it to the electrostatic simulation of low-dimensional semiconductor devices. The proposed approach is based on a fixed Cartesian grid and represents the geometry by level-set functions, avoiding mesh generation and remeshing even in the presence of interfaces. A snapping-back-to-grid strategy is used to control the small-cut-cell issue, while interface conditions are weakly enforced by a symmetric Nitsche formulation. The method is first validated on benchmark elliptic interface problems with different geometries and coefficient jumps, showing second-order accuracy for the solution and first-order accuracy for its gradient. As an application, we consider a graphene field-effect transistor described by a self-consistent drift-diffusion-Poisson model. The electrostatic potential is computed in a two-dimensional oxide/graphene/oxide structure, while charge transport in the graphene layer is modeled by one-dimensional bipolar drift-diffusion equations in the degenerate case and by including a field-dependent mobility model. The coupled nonlinear system is solved by a damped fixed-point iteration combined with a domain-decomposition treatment of the three-layer geometry. Numerical simulations reproduce transfer characteristics with a clear transition from an OFF state to an ON state and show the influence of the gate voltage on the two-dimensional electrostatic potential. The results also indicate that the effective thickness and discretization of the graphene layer affect the transverse potential profile, supporting the use of a full two-dimensional electrostatic description with explicit oxide/graphene interface conditions.

Computational solid mechanics relies on the discretisation of differential equations, requiring numerical differentiation and integration, often via quadrature. Recently, automatic differentiation (AD) has attracted interest as a tool to support these operations. This paper introduces the use of AD in nonlinear finite element analysis, focusing on stress updating, a key and repeatedly executed step in which stresses are computed from strain increments. Both implicit and explicit schemes require derivatives of constitutive models, which can be complex, error prone, and time consuming to derive analytically. We compare traditional analytical differentiation with automatic differentiation for implementing a hyperplastic Critical State model widely used in geotechnical engineering. Drained triaxial compression tests with varying overconsolidation ratios are simulated using both approaches. The results show that AD significantly simplifies the implementation of backward Euler stress integration by removing the need for manual derivation of derivatives. Despite this simplification, AD maintains high accuracy and robustness comparable to analytical approaches. The findings demonstrate that automatic differentiation can streamline the development of nonlinear material models, enabling efficient and reliable implementation of constitutive laws of arbitrary complexity. This opens the way for more flexible and maintainable computational mechanics codes.

Background: Over the last two decades, considerable effort in the applied mathematics community has been devoted to the development of polytopal methods, i.e., numerical methods for partial differential equations that support meshes more general than those used by traditional finite element methods. Despite their flexibility and attractive mathematical properties, polytopal methods have a steep implementation curve and, to date, no easily accessible open-source software framework has been available for their systematic implementation. Methods: In this article, we introduce poly_rust, a new software library implementing a domain-specific language (DSL) specifically designed for polytopal methods. poly_rust provides high-level constructs for expressing a wide range of such methods using a syntax that closely mirrors their mathematical formulation, thereby substantially lowering the implementation barrier for both newcomers and experienced users. Results: The library facilitates the rapid implementation, testing, and comparison of polytopal discretizations, making these methods more accessible to researchers, applied scientists, engineers, and students. In our vision, this will foster their dissemination beyond the numerical analysis community and encourage their adoption by practitioners in a broad range of application fields. It will also contribute to reproducibility of results, as DSL files are text files that can be easily made available, stored, and exchanged. Conclusions: poly_rust introduces, to our knowledge, the first DSL specifically tailored to polytopal methods. By combining mathematical expressiveness, ease of implementation, and an open-source framework, it provides a common platform for developing, sharing, reproducing, and comparing polytopal schemes, with the ambition of becoming a reference tool for the wider scientific community.

In this paper, we introduce a nonlocal nonlinear Schrödinger (NLS) model on periodic domains, and establish its well-posedness, conservation laws, local limit, dispersion properties, and develop a structure-preserving asymptotically compatible Fourier collocation method. We prove that for a fixed nonlocal horizon $\delta$, the model is globally well posed and conserves mass and nonlocal energy, and the nonlocal NLS solution converges to the local NLS solution with order $O(\delta^2)$ under suitable regularity assumptions. We also derive the dispersion relation and group velocity for plane waves, establish their $O(\delta^2)$ local limits, and characterize their behavior at high frequencies. The numerical method combines Crank–Nicolson time stepping with Fourier collocation and preserves mass and energy. The existence, uniqueness, and convergence of the numerical solutions are proved. In particular, the error of the nonlocal NLS solution is $O(\tau^2+N^{s-r})$ in the $H^s$ norm, uniformly with respect to the horizon. Moreover, its total $H^s$-error relative to the local NLS solution is $O(\delta^2+\tau^2+N^{s-r})$ without any coupling condition among $\delta$, $\tau$, and $N$, which proves asymptotic compatibility of the proposed method. Numerical experiments in one, two, and three dimensions are presented to verify the theoretical accuracy and discrete conservation, confirm convergence under independent variation of horizon and discretization parameters, and show how the horizon and kernel affect dispersive wave propagation.

A discrete entropy inequality is the principal nonlinear stability estimate available for systems of conservation laws, and evaluating it presupposes a physically admissible state. So far, however, the two have been secured separately. Entropy-stable schemes are almost always semi-discrete, are built around one selected entropy pair, and take for granted the positivity of density and pressure that makes the entropy well defined in the first place, whereas bound-preserving limiters keep the solution admissible but deliver no entropy estimate. For the special relativistic Euler equations, the two cannot be separated at all, since the conservative-to-primitive map is implicit, and an inadmissible state therefore has no entropy to correct. Here we construct high-order discontinuous Galerkin and finite volume schemes that, to our knowledge, for the first time, are entropy stable in the fully discrete sense for an arbitrary prescribed finite family of convex entropy pairs, a property we call multi-entropy stability, and are provably admissible wherever an entropy is evaluated. All of this is achieved by a single cellwise projection, and neither conservation nor the design order is lost. The construction rests on relativistic causality, which bounds every characteristic speed by the speed of light. Consequently, the numerical viscosity can be fixed once for all states and all equations of state, and one two-point building block then serves the whole entropy family. Since only the convexity of the admissible set and this speed bound are used, the same route remains open for related systems. Finally, in computations with four equations of state, the schemes retain high-order accuracy close to vacuum, produce no inadmissible state in strong shocks, near-vacuum shock–vortex interaction or jets with Lorentz factor above $70$, and confirm the monotone decay of every enforced discrete entropy.

We derive an explicit discrete energy identity for rational time discretizations generated by the first-subdiagonal Padé approximants of the exponential for solving linear seminegative problems. This work extends the diagonal Padé energy laws in [Z. Sun, Y. Wei, and K. Wu, SIAM J. Numer. Anal., 60 (2022)] to the first-subdiagonal family. The main new ingredient is an explicit Cholesky-type factorization of the energy coefficient matrix associated with the semi-inner-product terms in the discrete energy identity. The construction and proof of this factorization are nontrivial, since the matrix entries are alternating sums of Padé coefficients and the triangular factor has a parity-dependent factorial structure. We prove the factorization by reducing it to scalar rational identities and establishing them through finite product reductions and telescoping summations. Together with a \beta-coefficient cancellation, the factorization yields an exact discrete energy law that recovers the classical unconditional contractivity for linear seminegative problems. Numerical experiments adapted from the diagonal Padé energy-law setting illustrate the predicted order and verify the discrete dissipation identity.