Papers

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.

We study the effect of mixed precision on the numerical integration of the heat equation discretized with a Runge–Kutta method in time. A full space-time discretization is applied, which results in a very large and sparse linear system to be solved for the numerical approximations and the Runge–Kutta stages of all time steps. The linear system is solved by applying a suitable preconditioned iterative method that can be run in parallel. In order to speed up the solution process, the preconditioner is applied using a mixed precision framework. Sequential results show the robustness of the preconditioner, even when applied in mixed precision. Finally, we present numerical evidence of the improved performance of the mixed precision strategy when applied in a parallel environment, achieving up to a 50% reduction in CPU time.

Algebraic flux correction schemes for convection–diffusion–reaction equations rely on linearity preservation to avoid unnecessary artificial diffusion: the limiter is inactive whenever the discrete solution is affine. On quadrilateral and hexahedral meshes, this is insufficient, since the finite element space also contains bilinear functions that are reproduced exactly on axis-aligned meshes. We show that a mesh patch admits bilinearity preservation if and only if the central node lies in the convex hull of its neighbours after lifting each neighbour by the product of its coordinate differences, and we derive a sharp explicit value for the limiter parameter on such patches when the cell edges are parallel to the coordinate axes. We further characterize the meshes on which a mapped bilinear finite element space reproduces a bilinear function exactly, thereby establishing a fundamental limitation on what any limiter can achieve on general quadrilaterals. Numerical experiments confirm the theory: on axis-aligned grids, the proposed limiter reproduces a bilinear solution to the accuracy of the nonlinear iteration, whereas the linearity-preserving limiter does not. A sensitivity study shows that the parameter required for bilinearity preservation lies close to the upper end of the range for which the nonlinear problem can be solved by the fixed-point iteration considered here.

Penalizing a state constraint creates a singular localization problem: whole-space values may grow like the inverse penalty parameter, while numerical diffusion allows even an inward feedback to cross the boundary. For deterministic discounted optimal control, we show how the penalty itself supplies confinement that offsets this growth. With mesh size $h$ and penalty parameter $\varepsilon$, an inward barrier bounds the additional cost of numerical leakage by $O(h/\varepsilon)$ for a monotone centered-difference scheme with vanishing viscosity. Under $h\le\varepsilon$, an occupation estimate then yields $O(h)$ localization error with a sufficient box margin logarithmic in $1/h$ and independent of $\varepsilon$. We extend the result to bounded squared-distance penalties and combine it with residual-based evaluation bounds and discounted policy-error propagation. Under coupled refinement with vanishing penalization and discretization errors, explicit conditions on evaluation errors and greedy gaps ensure convergence of the neural value approximations to the constrained value, allowing measurable, nonunique greedy selectors. Reference calculations isolate leakage and localization across mesh and penalty scales, and distinguish evaluation, iteration, and approximation errors. An explicit cylindrical state-constraint solution provides a benchmark in arbitrary dimension, tested up to dimension twenty. Paired obstacle-navigation experiments illustrate why policy-value accuracy must be assessed alongside sampled residuals when comparing raw-residual and finite-grid-assisted neural evaluation.

Approximating the input-output behavior of a multivariable black-box function from limited data is challenging when blind to the importance of its inputs and their interactions. We introduce total sensitivity kernels (TSKs), a method based on families of weighted ANOVA kernels that learn and adapt to this multivariable structure. TSKs parameterize the weights on each multivariable component of the target function by factors for each input. We propose learning these factors directly from function evaluations by selecting the reproducing kernel Hilbert space (RKHS) in which the target function has minimum norm. Under suitable conditions, we show that this norm-minimization problem admits a unique solution, and we establish consistency of a finite-data formulation based on minimum-norm interpolation. The learned TSK factors characterize the participation of individual inputs across interactions and main effects, providing a kernel-dependent notion of input sensitivity related to total Sobol indices. Numerical experiments demonstrate that adapting the kernel to learned multivariable structure can substantially improve approximation accuracy over a standard product kernel.

Inverse problems are predominantly solved within the optimization-driven paradigm, which formulates the problem as objective-function minimization and approaches the solution by iterative search. Though universal, it suffers from limitations in efficiency, interpretability, and multi-parameter decoupling. This paper proposes a new paradigm—Structure-Driven Inversion (SDI). Instead of iterative search, SDI identifies and exploits the intrinsic structure of the problem to construct a solution method. It has two types of structure—mathematical structure and physical structure—and two corresponding driving inversion types: Mathematical Structure-Driven Inversion (MSDI) and Physical Structure-Driven Inversion (PSDI). The current representative method of MSDI is Mathematical Structure-Driven Pseudo-Inverse Inversion (MSDPII), which constructs a pseudo-inverse in the spectral domain via unitary diagonalization. The current representative methods of PSDI are Physical Structure-Driven Waveform Inversion Imaging (PSDWII), which performs inversion through virtual-source projection, and Physical Structure-Driven Back-Propagation (PSDBP), which implements three structural projections for deep neural networks via physical-system analogy. SDI is not a rejection but a complement and extension of Optimization-Driven Inversion (ODI). To the best of the author's knowledge, no existing study has systematically presented structure-driven inversion as an independent paradigm; this paper aims to fill that gap.

Reduced-order modeling of transport-dominated systems remains challenging because moving fronts and shocks are poorly represented by low-dimensional linear subspaces. We develop a structure-informed data-driven reduced-order model(ROM) for scalar hyperbolic conservation laws based on the kinetic defect formulation. This formulation separates the nonlinear dynamics into known characteristic transport and a kinetic entropy defect localized on the shock manifold. We exploit this structure by first removing the known transport from the solution snapshots. We then extract and register the remaining defect-driven dynamics in a shock-attached coordinate system. Separate ROMs are used to evolve the shock geometry and the registered defect-driven source. During prediction, the predicted shock geometry is used to inverse-register the learned defect-driven source, which advances the kinetic state and recovers the physical solution. Numerical examples in one and two spatial dimensions demonstrate accurate reconstruction and prediction of nonlinear transport with shocks, including evolution beyond the training interval, while accurately capturing the mass and entropy-dissipation behavior of the reference solution.

We propose an augmented high-order compact finite difference method for biharmonic equations with clamped boundary conditions and variable coefficients. Standard mixed-type formulations introduce an auxiliary variable, but its boundary values are unavailable, leaving the resulting discrete systems globally coupled and difficult to solve at large scales. Our key contribution is the development of a new augmented formulation that treats these unavailable boundary values as additional unknowns, reduces the global coupling to a lower-dimensional Schur complement system, and yields decoupled second-order subproblems. The Schur complement is solved by matrix-free GMRES, while the subproblems are handled by FFT-based fast solvers. The method achieves fourth-order accuracy using compact stencils, and has $O(n\log n)$ computational complexity, enabling the solution of the biharmonic equation with $1024^3$ degrees of freedom within several minutes. To the best of our knowledge, this level of computational efficiency has not previously been achieved in either the literature or practice. Using energy estimates and Fourier analysis, we derive a new $L^2$-estimate for Poisson equations with inexact Dirichlet boundary and then prove the convergence of the proposed scheme. We provide ample numerical experiments to confirm the accuracy, efficiency, and further apply the fast and accurate solver to triharmonic equations, high-wavenumber problems, Stokes flow, and plate bending problems.

We solve the four-dimensional Hamilton-Jacobi-Bellman (HJB) equation for two diffusively coupled Stuart-Landau-like oscillators to obtain full-state optimal feedback control. A sweep over coupling strength reveals a sharp change in the numerically optimal strategy: below a threshold coupling value, the controller steers the phase difference toward anti-phase while keeping both oscillators near the limit cycle, as reduced-order models would suggest. Above this threshold, the HJB solution changes qualitatively; the controller transiently collapses one oscillator's amplitude to near zero, thereby enabling large phase repositioning near the origin before rebuilding its amplitude. Direct gradient-based and stochastic optimization do not recover this lower-cost collapse trajectory from the initializations considered, suggesting that it occupies a region of the control landscape that is difficult to access by direct search. A joint sweep over nonisochronicity and coupling shows that collapse can occur even for an isochronous oscillator: phase repositioning near the origin can favor an off-cycle strategy. Nonisochronicity provides an additional energetic benefit through a phase-velocity surplus at small amplitude, quantitatively accounting for the observed reduction in control cost. Comparisons with uncoupled and coupled phase-reduced baselines show that phase models become increasingly inaccurate and cost significantly more energy for strong coupling. Results for coupled Van der Pol oscillators further demonstrate that exploitation of off-cycle dynamics is not specific to the Stuart-Landau-like oscillators.

Deep spline superposition networks face a tension between approximation order and stability across depth. We study approximation under a hard layerwise Lipschitz budget, and organise it around two quantities: the factorisation stability complexity of a given deep factorisation, and the budget-compatible approximation complexity of a discretisation operator. First, we solve exactly the finite-depth diagonal balancing problem for a fixed chain of nonnegative envelope matrices: the optimal uniform layer budget equals $\|M_{L-1}\cdots M_0\|_{\infty\to\infty}^{1/L}$, attained by an explicit one-pass minimiser, for rectangular layers, with a complete treatment of degeneracies and non-attainment. The optimum can be arbitrarily larger than the Lipschitz constant of the network itself, because passing to envelopes destroys sign cancellation. Second, we give a constructive spline discretisation theorem preserving the budget up to a controlled slack, with an explicit grid threshold. Conversely, for linear spline-valued operators that preserve the budget exactly, we prove budget-compatible minimax lower bounds on classes constrained simultaneously in the first and third derivative norms – a constraint pair that is forced by the problem and that rules out the usual scaling escapes. Finally, we show that the corresponding layer errors need not cancel under composition: for every operator of the class there is a stable depth-$L$ tower realising a constant fraction of the accumulated error, so the linear-in-depth accumulation of the upper bound is not a proof artefact.

Machine learning-based partial differential equations (PDEs) solvers have attracted significant attention in recent years. Most progress in this area has been driven by deep neural networks such as physics-informed neural networks (PINNs) and kernel method (such as physics-informed Gaussian Processes). We introduce a physics-informed random feature method for countering part of the spectral bias which PINN-based solvers are facing for a certain class of PDEs. Random feature method was originally proposed to approximate large-scale kernel machines and can be viewed as a specialized randomized neural network. Compared to other state-of-the-art PINN-based solvers which require a large number of collocation points, our proposed method reduces the computational complexity. In this paper, we develop a rigorous approximation error analysis and derive high-probability error bounds on the $H^1$ norm. We provide extensive numerical tests for verifying our theoretical guarantees on error decay rates, as well as several comparison tests to showcase our claimed capability for combating spectral bias in these deep learning based methods.

This work presents the first steps for simulating kinetic plasmas with collisions using the characteristic mapping method (CMM). The CMM is a semi-Lagrangian method that explores a semi-group structure to store diffeomorphic flow maps efficiently. Using the semi-group structure, individual submaps can be composed to relate the flow backward in time to its initial food point. The novelty of the presented work is handling the source term by storing sub-integrals that correspond to the individual sub-maps and allow efficient integration. Furthermore, we use the Lagrangian structure to avoid implicit time integration in the hydrodynamic regime, which is known to be stiff. We benchmark our method on the Boltzmann-BGK and Vlasov-Poisson-BGK equations and consider different test cases, the Sod shock tube problem and nonlinear Landau damping for different Knudsen numbers. We show third-order spatial and temporal convergence and illustrate the fine-scale zoom property of CMM for the bump-on-tail instability.

We introduce a new traffic flow model. The main idea is to describe acceleration effects in the Aw-Rascle-Zhang model through an additional driver-dependent equation, thereby allowing individual driver characteristics to evolve dynamically rather than being prescribed solely by the traffic density. The resulting system provides a hyperbolic description of traffic dynamics and admits an explicit analysis of its characteristic structure, Riemann problem, and entropy properties. To investigate the stability of the proposed model, we perform a Chapman-Enskog expansion up to second order and derive the corresponding diffusive and higher-order corrections. This analysis reveals how the stability properties depend on the underlying velocity and pressure functions. Several examples are presented to illustrate the theoretical results.

We study the temporal error growth of the Strang splitting method for the periodic cubic nonlinear Schrödinger (NLS) equation. The leading error is governed by a forced linearized equation, whose growth depends sharply on dimension and the sign of the nonlinearity. In the 1D defocusing case, we prove a uniform quadratic upper bound $C(1+T^2)\tau^2$ using the global Birkhoff transformation and the resulting degenerate structure of the linearized flow. In the higher-dimensional defocusing case, the exponential growth is constructed using arbitrarily small unstable standing waves. Moreover, to prove the higher-dimensional defocusing upper bound, we use the periodic Strichartz estimates of Killip and Vişan to show that $\int_0^T\|u(t)\|_{L^\infty}^2\,dt\le C(u_0)(1+T)$. This yields an exponential rate independent of the final time, despite possible growth of higher Sobolev norms. The exponential lower bound in the focusing case is obtained from the unstable linearized dynamics around a plane wave, with the Akhmediev breather providing the underlying mechanism. In addition, for 1D focusing case with small initial data, the appliction of the local Birkhoff transformation ensures that the error grows at most quadratically in time. Various numerical experiments confirm the sharp linear, quadratic, and exponential error growth rates.

Physically interpretable power take-off models support wave-energy development at every stage. For the seawater-pump class, few have been identified on a real machine and validated at the system level. This paper presents a block-modular grey-box model of a short-stroke seawater pump, identified from bench data alone. The pump decomposes into chamber, relief valve, piston seal and tip check valve, each structured by component physics, and the few parameters are identified block by block from terminal pressure, displacement, force and discharged-mass measurements, with a closed-cycle mass balance that decouples the degenerate seal and valve losses and closes to within 1% of the inflow. Over the sinusoidal campaign, the pressure-energised piston seal carries 84% of the lost mass and 97% of the lost energy, and the pump delivers about four-fifths of its input as high-pressure flow. Generalisation is demonstrated rather than assumed. With parameters fixed on ramps and sinusoids, the model reproduces 360 s sea-state records in free run to about 10% of range, captures 80% of several thousand peaks within 5 bar, transfers to a re-sprung valve on cracking strokes once its crack pressure is identified, follows thousand-cycle endurance records that drift below 0.1 bar, and reproduces the force and cycle-mean absorbed power, 0.1 to 5.7 kW, to a 4% median error. A practical-identifiability and uncertainty analysis confirms the terminal data resolve the retained parameters; thirty-seven repeated records set the experimental floor. The result is a compact, physically labelled model for design, energy assessment and model-based control of seawater-pump power take-offs.

The finite difference (FD) method is commonly used to approximate derivatives of smooth functions, with accuracy typically determined by stencil size and derivative order. However, certain centered stencils exhibit unexpectedly higher accuracy, a phenomenon known as superconvergence, which has been observed in practice but lacks rigorous explanation. We present a mathematical framework for superconvergence in centered FD approximations based on Taylor expansions of the truncation error and the resulting linear system for the FD coefficients. By analyzing symmetry properties of these coefficients and their interaction with the parity of the derivative order, we identify conditions under which higher-order error terms cancel. We show that superconvergence occurs for odd-order derivatives with even centered stencils and for even-order derivatives with odd centered stencils, while no superconvergence occurs for even derivatives with even centered stencils. Numerical experiments in MATLAB confirm the predicted convergence rates.

We present a new perspective on integral operator kernels given by a nonlinear function of the distance. This perspective is useful since it permits one to quasilinearize the nonlinear transformations acting on the tensor Haar expansion of the distance, d(x,y). This leads to a new representation comprised of a principal and residual term with desirable features; namely, enhanced local regularity of the principal component (and in certain situations rapid off-diagonal decay leading to O(N) storage) and a residual expansion of improved approximation as a consequence of enhanced global regularity. Numerical experiments on the potential kernel are included

Nested cubature rules, in which each refinement retains all previously used nodes and thus reuses earlier function evaluations, are natural in multilevel and adaptive integration. For every fixed $s>d/2$, we show that the equal-weight QMC integration rate on $\mathbb S^d$ is compatible with such nested point sets. In the subcritical range $d/2<s<d$, cumulative unions of geometrically growing QMC blocks yield nested QMC design sequences whose successive cardinality ratios converge to any prescribed $\rho>1$. At and above the critical index $s=d$, where the block-averaging estimate no longer yields the optimal rate, we prove an equal-weight completion theorem based on low-frequency discrepancy cancellation. A block-sensitive estimate sharpens the iteration and yields nested QMC designs for all $s\ge d$ with $N_{j+1}\lesssim(j+1)^{2s/d-1}N_j$. Every prescribed finite point set also admits an optimal-rate completion.

In this work, we propose an unfitted finite element scheme to approximate the solution of the heat equation on moving domains. We use the {\phi}-FEM paradigm in which the computational domain is described implicitly by a level-set function {\phi}. This function is incorporated into the variational formulation in order to enforce the boundary conditions, which allows the use of unfitted meshes in space. Such a strategy avoids the need for remeshing at each time step and makes it possible to handle complex geometrical evolutions. Moreover, the {\phi}-FEM approach has the advantage of being simple to implement within standard finite element libraries. We introduce a fully discrete scheme that combines the {\phi}-FEM spatial discretization with the lowest-order discontinuous Galerkin method in time. Under regularity assumptions on the level-set function and a mild restriction linking the time step to the mesh size, we establish an optimal a priori error estimate in the L2(0,T;H1) norm. Finally, we present several numerical experiments that confirm this convergence rate, exhibit a second-order convergence in the L{\infty}(0,T;L2) norm, and illustrate the robustness and accuracy of the proposed method.

Asynchronous iterative methods are attractive for large-scale parallel computing because they reduce synchronization and communication overhead. Existing convergence-rate analyses, however, have primarily focused on shared memory implementations, whereas distributed memory systems introduce more general and potentially inconsistent communication delays. In this work, we revisit asynchronous Jacobi and randomized Gauss–Seidel (RGS) methods for symmetric positive definite linear systems from a unified perspective. We first introduce a general asynchronous model that encompasses both shared memory and distributed memory settings and expresses the two methods through a common coordinate-update framework. We then establish linear convergence in expectation under an explicit stability condition. The resulting convergence bound depends algebraically on the delay through the quantity $\sqrt{\rho\tau}+\rho\tau$, where $\tau$ is the maximum communication delay and $\rho$ reflects the communication pattern of the underlying parallel implementation. In particular, under an appropriate scaling regime with $\rho\tau=O(1)$, the guaranteed per-iteration convergence rate has the same asymptotic order as that of synchronous RGS. These results provide a unified framework for quantifying the effect of asynchronicity on asynchronous Jacobi/RGS methods in both shared and distributed memory environments.

The fast Fourier transform (FFT) and its real counterparts, the fast algorithms of the discrete cosine transform (DCT) and the discrete sine transform (DST), are frequently employed for the numerical evaluation of the Fourier cosine and Fourier sine transform, thereby acting as parametric quadrature rules based on finitely many samples of a given continuous function. But how accurate are the obtained results? In this paper we study the error occurring if the DCT of type I (DCT-I) and the DST of type I (DST-I) are applied to compute $\int_0^{\infty} f(x)\cos(2\pi xv)\d x$ and $\int_0^{\infty} f(x)\sin(2\pi xv)\d x$ for $f \in L^1(\Rp) \cap C(\Rp)$ and $v \in \Rp$, as well as the error occurring if the Chebyshev coefficients $\frac{2}{\pi}\int_0^{\pi} h(\cos\theta)\,\cos(n\theta)\d\theta$, $n \in \Np$, of a function $h \in C(I)$ on $I := [-1,1]$ are computed by the DCT-I. We present practicable, explicit error estimates under polynomial, exponential, and mixed decay conditions. For the Fourier cosine and sine transform, rational functions lead to a geometric decay in the frequency parameter $P$ combined with an algebraic plateau of order $L^{-3/2}$ in the sampling parameter $L$, while for meromorphic functions such as $\operatorname{sech}(\pi\,\cdot)$ the critical geometric rate is attained in both parameters. For Chebyshev approximation of a rational function with simple poles in $\C \setminus I$, the quadrature error has the critical geometric rate $\rho_0^{-P}$ with an explicit constant, where $\rho_0 > 1$ is the parameter of the largest Bernstein ellipse of analyticity. Several numerical examples illustrate the theory.

For every reproducing kernel Hilbert space $\mathcal{H}_K$ with bounded kernel $K$, we prove that the linear sampling widths $g_m^{\text{lin}}$ and Gelfand widths $c_n$ in the uniform norm satisfy $$ g_m^{\text{lin}}(B_{\mathcal H_K})_\infty \leq \frac{m+1}{m-n+1}\, c_n(B_{\mathcal H_K})_\infty\quad , \quad m\ge n. $$ In a certain sense, the leading factor $(m+1)/(m-n+1)$ on the right-hand side is sharp. We transfer the problem into the selection of a maximal volume subset of kernel translates (kernel $D$-optimal design). Apart from basic linear algebra, our proof employs low rank approximation techniques and a version of the Eckart-Young-Mirsky theorem. The same argument gives corresponding bounds for the selection of a reduced basis in Hilbert spaces.

We consider partial differential equation (PDE) inverse problems from sparse, noisy, and corrupted observations, with the aim of recovering unknown coefficient fields and associated state variables in a mesh-free setting. Under such sparse and corrupted observations, standard physics-informed and generative approaches typically treat all samples indiscriminately and therefore lack a principled mechanism for reconciling physical laws with contaminated data. We address this difficulty with a two-stage flow-matching framework. In the first stage, we develop physics-guided conditional flow matching (PG-CFM), which incorporates strong-form PDE information through residual regularization along the generative trajectories together with intermittent global collocation constraints. In the second stage, we introduce energy-regularized flow matching (ERFM), which fine-tunes the Stage-1 model by assigning each observation a physics–data energy score from a frozen teacher and reweighting the flow-matching objective to reduce the influence of high-energy, PDE-inconsistent samples. We show that the resulting Stage-2 objective is equivalent to flow matching under a teacher-induced reweighted data distribution, which gives a population-level interpretation of the robustness mechanism. Numerical experiments on several inverse benchmarks, including Poisson and Navier–Stokes problems, show that the proposed framework yields more accurate coefficient recovery than robust PINN variants and competing generative baselines.

We propose a mixed finite element discretisation for the incompressible Navier-Stokes equations that preserves the evolution laws of both energy and enstrophy, in a stronger sense than previous discretisations. In two dimensions, the evolution law for enstrophy only permits dissipation for thermodynamically isolated systems, leading to a Reynolds-number-independent bound on the velocity gradient that naturally stabilises the scheme, even on severely under-resolved meshes. In three dimensions, the scheme preserves both dissipation and the generation of enstrophy through vortex stretching. We enforce these evolution laws by systematically introducing auxiliary variables into the discretisation. While conforming implementations of these schemes require discrete Stokes complexes with enhanced regularity, we introduce both (i) equivalent reparametrisations and (ii) penalty formulations that require only the typical curl- and div-conforming spaces from the standard discrete de Rham complex. The scheme handles different types of boundary conditions and curved domains. The robust stabilisation properties of the proposed scheme are demonstrated through numerical simulations of a shear flow, a spherical vortex, and flow past an obstacle. We observe numerically that preserving the discrete evolution of enstrophy in this way has a strong stabilising effect on the numerical solution, especially in two dimensions.

We introduce a fully–discrete variational approximation of Mather measures and sets for Tonelli Lagrangians on the flat torus, together with a numerical procedure for approximating the entire Mather set. The scheme is based on a fully–discrete Lax–Oleinik operator with integer winding labels. We prove an $O(\tau+h/\tau)$ error estimate for the critical value, convergence of critical solutions, and a finite-dimensional characterization of fully–discrete Mather measures. Accumulation points of the reconstructed minimizing measures are continuous Mather measures, while the supports satisfy complementary upper and lower convergence results involving the Mañé and Mather sets. To avoid the selection of only some minimizing components by exact discrete minimizers, we introduce a mass-threshold approximation based on almost-minimizing holonomic measures. This yields a finite-dimensional constrained optimization procedure designed to recover the whole Mather set.

We propose a dissipative, irreversible damage model for anisotropic media in large deformations to address the damage process of the myocardium following a cardiac infarction. We model damage to the collagen in the cleavage planes resulting from an increased load to the passive tissue. Starting from variational principles, we derive a quasi-static differential model governing the irreversible evolution of the damage. We apply our model to two test cases. First, a rectangular passive slab geometry to which an indenter-like load is applied. Then, we couple our model with a model for the numerical simulation of left ventricular electromechanics. Our results show that a redistribution of the load within a left ventricle occurring following a cardiac infarction can produce the damage to the collagen scaffold in the infarcted area consistent with experimental results. This is indicative of the fact that passive load-dependent dissipative phenomena are implicated in post-infarction remodelling.

In this paper we propose augmented Lagrangian finite element methods for cavitation in the Reynolds and Stokes models of lubrication, for which the discrete complementarity conditions hold exactly and the augmentation parameter drops out of the method. For the Reynolds equation this follows from nodal quadrature in a mixed piecewise linear method, for which we prove the classical first-order error estimate. We also determine the computable stability threshold of a multiplier-free stabilised alternative. For Stokes flow we show that the constrained scalar is the mechanical pressure when the deviatoric stress, with zero bulk viscosity, is used, but not with the customary incompressible stress. We discretise with a jump-stabilised Crouzeix-Raviart element and piecewise constant pressure, and prove well-posedness, stability in two and three dimensions, and a first-order error estimate. We present numerical examples which verify these properties and show that close pressure profiles in the two models do not imply close cavity predictions, the latter also being sensitive to the end conditions of the computational domain.

We propose a hybridized staggered discontinuous Galerkin–mixed finite element method for the strain gradient elasticity model, a fourth-order singularly perturbed problem governed by a material length scale parameter $\iota$ and the Lamé constants $\lambda$ and $\mu$. Introducing the total stress together with the scaled Cauchy and hyper stresses as auxiliary unknowns, we recast the model as a first-order system and discretize the displacement and the total stress by staggered discontinuous Galerkin spaces, the two scaled stresses by Raviart–Thomas pairs, with symmetry imposed strongly on all three stresses. The higher-order Dirichlet boundary condition enters the variational formulation naturally, so that the scheme avoids the numerical boundary layer that limits displacement-based methods. By establishing an inf-sup condition on the symmetric subspace and a discrete Korn inequality, we prove algebraic stability and optimal convergence, both uniform in $\iota$ and $\lambda$. We further develop a hybridized scheme in which local static condensation leaves only two multipliers on the mesh skeleton as globally coupled unknowns. Numerical experiments confirm the predicted convergence rates and the parameter robustness.

We present shape calculus techniques and a scalable automatic shape differentiation framework for multi-phase topology optimisation on unfitted discretisations defined by several level-set functions. First we establish general, exact shape calculus expressions in the discrete case for multi-phase systems by leveraging concepts from convex geometry. To complement this theoretical foundation, we introduce an open-source multi-phase automatic shape differentiation framework based on polytopal cutting. This computational framework is validated against both finite differences and our established exact expressions, matching the latter to near machine precision. Furthermore, the proposed automatic shape differentiation is scalable across distributed computing environments, demonstrating near ideal weak scaling up to 1.65 billion finite elements across 13,824 computer cores. We demonstrate our implementation by solving unfitted multi-phase topology optimisation problems for anisotropic diffusion, linear elasticity, and fluid–structure interaction. Together, these theoretical and computational contributions provide a robust and accessible foundation for advancing multi-phase topology optimisation using unfitted finite element methods. In particular, the methods enable the solution of topology optimisation problems involving multi-phase and multi-physics systems with non-trivial boundary conditions. The open-source software is available at https://github.com/zjwegert/GridapTopOpt.jl.

First-arrival traveltime tomography provides a computationally efficient means of imaging subsurface velocity structure, while boundary-conforming deformed grids allow rugged topography to be represented without abandoning a logically structured mesh. In matrix-free Gauss-Newton inversion, however, the linearized full-metric Eikonal transport and its transpose must be applied repeatedly at a fixed background model. Directional sweeping consequently revisits the same state-dependent dependency structure for every Krylov right-hand side. We recast the frozen tangent transport as a directed graph and distinguish loss of the background-traveltime ordering from genuine algebraic cyclicity. Strongly connected components identify the irreducible part of the transport, whereas the remaining dependencies admit exact scalar substitution after reordering. The condensation graph and local block factorizations are constructed once for each frozen source state and then reused for both primal and transpose applications. Numerical experiments on deformed grids and a three-dimensional tomography problem show that the resulting block-triangular formulation preserves the frozen sensitivity action while substantially reducing the repeated transport cost. By eliminating this dominant inner-iteration cost, the proposed method markedly accelerates Gauss-Newton computation and makes full-metric matrix-free inversion on deformed grids practical at much lower cost.