Preserving the admissible set of the ideal magnetohydrodynamics (MHD) equations is important not only for producing physically meaningful numerical solutions, but more importantly for achieving robust computations. In this paper, we develop an optimization-based limiter to enforce admissibility while preserving global conservation and accuracy. For an easy and efficient projection, we decompose the admissible set into slices parameterized by the magnetic energy, so that the MHD projection reduces to a one-dimensional minimization, which can be solved efficiently by the Brent method. The splitting method can be used to efficiently solve the global minimization problem of the optimization-based limiter, which can be used to enforce cell average admissibility in discontinuous Galerkin (DG) schemes, and pointwise admissibility can be further enforced by the Zhang-Shu positivity-preserving limiter. We apply the limiter to high-order DG schemes and present numerical results for a few representative MHD problems.

We extend the tensor-product direct solver from the Laplacian to the Schrödinger operator -Δ+ V. When the potential V_1 is separable, the operator -Δ+ V_1 is inverted or exponentiated at cost O(N^1+1/d) in d dimensions via per-axis eigendecomposition. On a single NVIDIA A100 GPU, this costs less than one second for 10^9 degrees of freedom in 3D. For non-separable potentials V = V_1 + V_2, the same solver provides a preconditioner (-Δ+ V_1)^-1 for the preconditioned conjugate gradient (PCG) method and a propagator for operator-splitting time integrators. For bounded V_2, we prove that the preconditioned operator has a bounded condition number and a clustered spectrum with at most finitely many outlier eigenvalues, independently of the mesh size, and also independently of the domain size when V_1 is a confining potential. This explains the mesh- and domain-independent PCG iteration counts observed in practice. We apply this method to ground state computation via inverse iteration for linear problems and via the a_u gradient flow for Gross–Pitaevskii energy in 3D, and also Hamiltonian simulation via the approximated qHOP and Magnus-2 splitting methods from 3D to 9D on a single NVIDIA GH200 GPU.

In low-rank numerical methods for Vlasov dynamics, the SVD-type truncation procedure may introduce negative entries into the numerical solution. Such negative values are unphysical because the solution is a probability distribution function. We design optimization-based post-processing algorithms to recover nonnegativity while preserving the macroscopic quantities (density, momentum, and energy) pointwise. The preservation of the macroscopic quantities is written as an orthogonality constraint on the correction term. For a convex formulation based on squared nuclear norm minimization, we show that the proximal operator with the orthogonality constraint is characterized by an implicit singular value thresholding equation, and the threshold can be computed efficiently by bisection. Based on this result, we develop five algorithms for the convex formulation: Douglas–Rachford splitting, restarted dual FISTA, restarted dual accelerated gradient descent, dual PR+ conjugate gradient, and dual L-BFGS. We also consider a non-convex formulation with an explicit rank constraint and develop a tangent-space accelerated alternating projection algorithm that only requires a 2r × 2r SVD per iteration. Numerical results for a Landau damping test case show that the proposed algorithms give comparable correction quality. Among them, the tangent-space accelerated alternating projection is the most cost-efficient, increasingly so as the problem size grows. We further demonstrate the correction as a positivity limiter inside a time-dependent conservative low-rank Vlasov solver, where it removes the negativity introduced by the SVD-type truncation while preserving the conserved mass, momentum, and energy.

High-order numerical methods for conservation laws can generate spurious oscillations near discontinuities. We propose to suppress these oscillations by post-processing the numerical solution with a total variation (TV) denoising step that reduces the total variation of the nodal values. For a general analysis operator, the resulting discrete minimization problem need not satisfy the submodularity property that existing exact max-flow algorithms require, so we instead compute the TV minimizer exactly using a differential inclusion algorithm. In one dimension, the differential inclusion algorithm computes the TV minimizer in finitely many steps, requires no tuning of algorithmic parameters, and preserves the total mass. In two dimensions, we apply the 1D algorithm dimension by dimension. We test the method as a post-processing limiter for Fourier pseudospectral methods and a fifth-order finite difference scheme applied to scalar conservation laws, compressible Euler equations, and the two-dimensional incompressible Euler equations.

We introduce effective splitting methods for implementing optimization-based limiters to enforce the invariant domain in gas dynamics in high order accurate numerical schemes. The key ingredients include an easy and efficient explicit formulation of the projection onto the invariant domain set, and also proper applications of the classical Douglas-Rachford splitting and its more recent extension Davis-Yin splitting. Such an optimization-based approach can be applied to many numerical schemes to construct high order accurate, globally conservative, and invariant-domain-preserving schemes for compressible flow equations. As a demonstration, we apply it to high order discontinuous Galerkin schemes and test it on demanding benchmarks to validate the robustness and performance of both ℓ^1-norm minimization limiter and ℓ^2-norm minimization limiter.

We compute the ground state u of the Gross–Pitaevskii equation (GPE) via Wasserstein gradient descent in diffeomorphism space. We represent the density ρ=u^2 as the push-forward of a fixed reference measure through a parameterized transport map T_θ, realized by a boundary-preserving Neural ODE. The Wasserstein gradient flow on probability densities then lifts to natural gradient descent in the finite-dimensional parameter space, with metric tensor given by the pullback of the Wasserstein metric. The method is entirely mesh-free and preserves the unit-mass constraint without normalization. We present numerical experiments in dimensions d=1,2,3 and demonstrate that the parameterized Wasserstein gradient flow (PWGF) output can be used to initialize the H^1 Sobolev gradient flow, reducing the initial energy gap by a factor of 7 in 2D and 4.5 in 3D compared to trivial initial conditions.

In this paper, we develop bound-preserving techniques for the Runge–Kutta (RK) discontinuous Galerkin (DG) method with compact stencils (cRKDG method) for hyperbolic conservation laws. The cRKDG method was recently introduced in [Q. Chen, Z. Sun, and Y. Xing, SIAM J. Sci. Comput., 46: A1327–A1351, 2024]. It enhances the compactness of the standard RKDG method, resulting in reduced data communication, simplified boundary treatments, and improved suitability for local time marching. This work improves the robustness of the cRKDG method by enforcing desirable physical bounds while preserving its compactness, local conservation, and high-order accuracy. Our method is extended from the seminal work of [X. Zhang and C.-W. Shu, J. Comput. Phys., 229: 3091–3120, 2010]. We prove that the cell average of the cRKDG method at each RK stage preserves the physical bounds by expressing it as a convex combination of three types of forward-Euler solutions. A scaling limiter is then applied after each RK stage to enforce pointwise bounds. Additionally, we explore RK methods with less restrictive time step sizes. Because the cRKDG method does not rely on strong-stability-preserving RK time discretization, it avoids its order barriers, allowing us to construct a four-stage, fourth-order bound-preserving cRKDG method. Numerical tests on challenging benchmarks are provided to demonstrate the performance of the proposed method.

Admissible states in hyperbolic systems and related equations often form a convex invariant domain. Numerical violations of this domain can lead to loss of hyperbolicity, resulting in illposedness and severe numerical instabilities. It is therefore crucial for numerical schemes to preserve the invariant domain to ensure both physically meaningful solutions and robust computations. For complex systems, constructing invariant-domain-preserving (IDP) schemes is highly nontrivial and particularly challenging for high-order accurate methods. This paper presents a comprehensive survey of IDP schemes for hyperbolic and related systems, with a focus on the most popular approaches for constructing provable IDP schemes. We first give a systematic review of the fundamental approaches for establishing the IDP property in first-order accurate schemes, covering finite difference, finite volume, finite element, and residual distribution methods. Then we focus on two widely used and actively developed classes of high order IDP schemes as well as their recent developments, most of which have emerged in the past decade. The first class of methods seeks an intrinsic weak IDP property in high-order schemes and then designs polynomial limiters to enforce a strong IDP property at the points of interest. This generic approach applies to high-order finite volume and discontinuousGalerkin schemes. The second class is based on the flux limiting approaches, which originated from the flux-corrected transport method and can be adapted to a broader range of spatial discretizations, including finite difference and continuous finite element methods. In this survey, we elucidate the main ideas in the construction of IDP schemes, provide some new perspectives and insights, with extensive examples, and numerical experiments in gas dynamics and magnetohydrodynamics.

We prove an explicit local linear rate for ADMM solving the isotropic Total Variation (TV) norm compressed sensing problem in multiple dimensions, by analyzing the auxiliary variable in the equivalent Douglas-Rachford splitting on a dual problem. Numerical verification on large 3D problems and real MRI data will be shown. Though the proven rate is not sharp, it is close to the observed ones in numerical tests.

For smooth optimization problems with a Hermitian positive semidefinite fixed-rank constraint, we consider three existing approaches including the simple Burer–Monteiro method, and Riemannian optimization over quotient geometry and the embedded geometry. These three methods can be all represented via quotient geometry with three Riemannian metrics g^i(· , · ) (i=1,2,3) . By taking the nonlinear conjugate gradient method (CG) as an example, we show that CG in the factor-based Burer–Monteiro approach is equivalent to Riemannian CG on the quotient geometry with the Bures–Wasserstein metric g^1 . Riemannian CG on the quotient geometry with the metric g^3 is equivalent to Riemannian CG on the embedded geometry. For comparing the three approaches, we analyze the condition number of the Riemannian Hessian near a minimizer under the three different metrics. Under certain assumptions, the condition number from the Bures–Wasserstein metric g^1 is significantly worse than the other two metrics. Numerical experiments show that the Burer–Monteiro CG method has obviously slower asymptotic convergence rate either when the minimizer has a large condition number or when it is rank deficient, which is consistent with the condition number analysis.

We first briefly review some recently proven new results about Q^k spectral element method for second order linear PDEs, including its order of accuracy as a finite difference method in ℓ ^2 -norm and monotonicity, both of which are special properties of Q^k spectral element method on structured meshes. We discuss some extensions or applications of these two special properties, including the accuracy for the Helmholtz equation and applications of monotone discrete Laplacian to a semi-linear problem. In particular, the Q^2 spectral element method gives a fourth order accurate monotone discrete Laplacian, with which one can obtain explicit convergence rates of Picard and Newton iterations for solving a special second order semilinear PDE.

This paper studies the numerical approximation of the ground state of the Gross-Pitaevskii (GP) eigenvalue problem with a fully discretized Sobolev gradient flow induced by the H1 norm. For the spatial discretization, we consider the finite element method with quadrature using Pk basis on a simplicial mesh and Qk basis on a rectangular mesh. We prove the global convergence to a critical point of the discrete GP energy, and establish a local exponential convergence to the ground state under the assumption that the linearized discrete Schrodinger operator has a positive spectral gap. We also show that for the P1 finite element discretization with quadrature on an unstructured shape regular simplicial mesh, the eigengap satisfies a mesh-independent lower bound, which implies a mesh-independent local convergence rate for the proposed discrete gradient flow. Numerical experiments with discretization by high-order Qk spectral element methods in two and three dimensions are provided to validate the efficiency of the proposed method.

For high-order accurate schemes such as discontinuous Galerkin (DG) methods solving Fokker–Planck equations, it is desired to efficiently enforce positivity without losing conservation and high-order accuracy, especially for implicit time discretizations. We consider an optimization-based positivity-preserving limiter for enforcing positivity of cell averages of DG solutions in a semi-implicit time discretization scheme, so that the point values can be easily enforced to be positive by a simple scaling limiter on the DG polynomial in each cell. The optimization can be efficiently solved by a first-order splitting method with nearly optimal parameters, which has an O ( N ) computational complexity and is flexible for parallel computation. Numerical tests are shown on some representative examples to demonstrate the performance of the proposed method.

In this paper, we introduce a quasi-Newton method optimized for efficiently solving quasi-linear elliptic equations and systems, with a specific focus on GPU-based computation. By approximating the Jacobian matrix with a combination of linear Laplacian and simplified nonlinear terms, our method reduces the computational overhead typical of traditional Newton methods while handling the large, sparse matrices generated from discretized PDEs. We also provide a convergence analysis demonstrating local convergence to the exact solution under optimal choices for the regularization parameter, ensuring stability and efficiency in each iteration. Numerical experiments in two- and three-dimensional domains validate the proposed method’s robustness and computational gains with tensor product implementation. This approach offers a promising pathway for accelerating quasi-linear elliptic equations and systems solvers, expanding the feasibility of complex simulations in physics, engineering, and other fields leveraging advanced hardware capabilities.

Polynomial approximations of functions are widely used in scientific computing. In certain applications, it is often desired to require the polynomial approximation to be non-negative (resp. non-positive), or bounded within a given range, due to constraints posed by the underlying physical problems. Efficient numerical methods are thus needed to enforce such conditions. In this paper, we discuss effective numerical algorithms for polynomial approximation under non-negativity constraints. We first formulate the constrained optimization problem, its primal and dual forms, and then discuss efficient first-order convex optimization methods, with a particular focus on high dimensional problems. Numerical examples are provided, for up to 200 dimensions, to demonstrate the effectiveness and scalability of the methods.

In many applications, it is desired to obtain extreme eigenvalues and eigenvectors of large Hermitian matrices by efficient and compact algorithms. In particular, orthogonalization-free methods are preferred for large-scale problems for finding eigenspaces of extreme eigenvalues without explicitly computing orthogonal vectors in each iteration. For the top $p$ eigenvalues, the simplest orthogonalization-free method is to find the best rank-$p$ approximation to a positive semi-definite Hermitian matrix by algorithms solving the unconstrained Burer-Monteiro formulation. We show that the nonlinear conjugate gradient method for the unconstrained Burer-Monteiro formulation is equivalent to a Riemannian conjugate gradient method on a quotient manifold with the Bures-Wasserstein metric, thus its global convergence to a stationary point can be proven. Numerical tests suggest that it is efficient for computing the largest $k$ eigenvalues for large-scale matrices if the largest $k$ eigenvalues are nearly distributed uniformly.