We propose two Hybrid High-Order (HHO) methods for the incompressible Navier-Stokes equations and investigate their robustness with respect to the Reynolds number. While both methods rely on a HHO formulation of the viscous term, the pressure-velocity coupling is fundamentally different, up to the point that the two approaches can be considered antithetical. The first method is kinetic energy preserving, meaning that the skew-symmetric discretization of the convective term is guaranteed not to alter the kinetic energy balance. The approximated velocity fields exactly satisfy the divergence free constraint and continuity of the normal component of the velocity is weakly enforced on the mesh skeleton, leading to H-div conformity. The second scheme relies on Godunov fluxes for pressure-velocity coupling: a Harten, Lax and van Leer (HLL) approximated Riemann Solver designed for cell centered formulations is adapted to hybrid face centered formulations. The resulting numerical scheme is robust up to the inviscid limit, meaning that it can be applied for seeking approximate solutions of the incompressible Euler equations. The schemes are numerically validated performing steady and unsteady two dimensional test cases and evaluating the convergence rates on h-refined mesh sequences. In addition to standard benchmark flow problems, specifically conceived test cases are conducted for studying the error behaviour when approaching the inviscid limit.
Higher-order time integration methods that unconditionally preserve the positivity and linear invariants of the underlying differential equation system cannot belong to the class of general linear methods. This poses a major challenge for the stability analysis of such methods since the new iterate depends nonlinearly on the current iterate. Moreover, for linear systems, the existence of linear invariants is always associated with zero eigenvalues, so that steady states of the continuous problem become non-hyperbolic fixed points of the numerical time integration scheme. Altogether, the stability analysis of such methods requires the investigation of non-hyperbolic fixed points for general nonlinear iterations. Based on the center manifold theory for maps we present a theorem for the analysis of the stability of non-hyperbolic fixed points of time integration schemes applied to problems whose steady states form a subspace. This theorem provides sufficient conditions for both the stability of the method and the local convergence of the iterates to the steady state of the underlying initial value problem. This theorem is then used to prove the unconditional stability of the MPRK22($\alpha$)-family of modified Patankar-Runge-Kutta schemes when applied to arbitrary positive and conservative linear systems of differential equations. The theoretical results are confirmed by numerical experiments.
We develop structure-preserving finite element methods for the incompressible, resistive Hall magnetohydrodynamics (MHD) equations. These equations incorporate the Hall current term in Ohm's law and provide a more appropriate description of fully ionized plasmas than the standard MHD equations on length scales close to or smaller than the ion skin depth. We introduce a stationary discrete variational formulation of Hall MHD that enforces the magnetic Gauss's law exactly (up to solver tolerances) and prove the well-posedness and convergence of a Picard linearization. For the transient problem, we present time discretizations that preserve the energy and magnetic and hybrid helicity precisely in the ideal limit for two types of boundary conditions. Additionally, we present an augmented Lagrangian preconditioning technique for both the stationary and transient cases. We confirm our findings with several numerical experiments.
Parameter reconstructions are indispensable in metrology. Here, on wants to explain $K$ experimental measurements by fitting to them a parameterized model of the measurement process. The model parameters are regularly determined by least-square methods, i.e., by minimizing the sum of the squared residuals between the $K$ model predictions and the $K$ experimental observations, $\chi^2$. The model functions often involve computationally demanding numerical simulations. Bayesian optimization methods are specifically suited for minimizing expensive model functions. However, in contrast to least-square methods such as the Levenberg-Marquardt algorithm, they only take the value of $\chi^2$ into account, and neglect the $K$ individual model outputs. We introduce a Bayesian target-vector optimization scheme that considers all $K$ contributions of the model function and that is specifically suited for parameter reconstruction problems which are often based on hundreds of observations. Its performance is compared to established methods for an optical metrology reconstruction problem and two synthetic least-squares problems. The proposed method outperforms established optimization methods. It also enables to determine accurate uncertainty estimates with very few observations of the actual model function by using Markov chain Monte Carlo sampling on a trained surrogate model.
The canonical polyadic decomposition (CPD) is a fundamental tensor decomposition which expresses a tensor as a sum of rank one tensors. In stark contrast to the matrix case, with light assumptions, the CPD of a low rank tensor is (essentially) unique. The essential uniqueness of CPD makes this decomposition a powerful tool in many applications as it allows for extraction of component information from a signal of interest. One popular algorithm for algebraic computation of a CPD is the generalized eigenvalue decomposition (GEVD) which selects a matrix subpencil of a tensor, then computes the generalized eigenvectors of the pencil. In this article, we present a simplification of GEVD which improves the accuracy of the algorithm. Surprisingly, the generalized eigenvector computation in GEVD is in fact unnecessary and can be replaced by a QZ decomposition which factors a pair of matrices as a product of unitary and upper triangular matrices. Computing a QZ decomposition is a standard first step when computing generalized eigenvectors, so our algorithm can been seen as a direct simplification of GEVD.
We propose, study, and compute solutions to a class of optimal control problems for hyperbolic systems of conservation laws and their viscous regularization. We take barotropic compressible Navier--Stokes equations (BNS) as a canonical example. We first apply the entropy--entropy flux--metric condition for BNS. We select an entropy function and rewrite BNS to a summation of flux and metric gradient of entropy. We then develop a metric variational problem for BNS, whose critical points form a primal-dual BNS system. We design a finite difference scheme for the variational system. The numerical approximations of conservation laws are implicit in time. We solve the variational problem with an algorithm inspired by the primal-dual hybrid gradient method. This includes a new method for solving implicit time approximations for conservation laws, which seems to be unconditionally stable. Several numerical examples are presented to demonstrate the effectiveness of the proposed algorithm.
The analysis of structure-preserving numerical methods for the Poisson--Nernst--Planck (PNP) system has attracted growing interests in recent years. In this work, we provide an optimal rate convergence analysis and error estimate for finite difference schemes based on the Slotboom reformulation. Different options of mobility average at the staggered mesh points are considered in the finite-difference spatial discretization, such as the harmonic mean, geometric mean, arithmetic mean, and entropic mean. A semi-implicit temporal discretization is applied, which in turn results in a non-constant coefficient, positive-definite linear system at each time step. A higher order asymptotic expansion is applied in the consistency analysis, and such a higher order consistency estimate is necessary to control the discrete maximum norm of the concentration variables. In convergence estimate, the harmonic mean for the mobility average, which turns out to bring lots of convenience in the theoretical analysis, is taken for simplicity, while other options of mobility average would also lead to the desired error estimate, with more technical details involved. As a result, an optimal rate convergence analysis on concentrations, electric potential, and ionic fluxes is derived, which is the first such results for the structure-preserving numerical schemes based on the Slotboom reformulation. It is remarked that the convergence analysis leads to a theoretical justification of the conditional energy dissipation analysis, which relies on the maximum norm bounds of the concentration and the gradient of the electric potential. Some numerical results are also presented to demonstrate the accuracy and structure-preserving performance of the associated schemes.
Group synchronization asks to recover group elements from their pairwise measurements. It has found numerous applications across various scientific disciplines. In this work, we focus on orthogonal and permutation group synchronization which are widely used in computer vision such as object matching and structure from motion. Among many available approaches, the spectral methods have enjoyed great popularity due to their efficiency and convenience. We will study the performance guarantees of the spectral methods in solving these two synchronization problems by investigating how well the computed eigenvectors approximate each group element individually. We establish our theory by applying the recent popular~\emph{leave-one-out} technique and derive a~\emph{block-wise} performance bound for the recovery of each group element via eigenvectors. In particular, for orthogonal group synchronization, we obtain a near-optimal performance bound for the group recovery in presence of additive Gaussian noise. For permutation group synchronization under random corruption, we show that the widely-used two-step procedure (spectral method plus rounding) can recover all the group elements exactly if the SNR (signal-to-noise ratio) is close to the information theoretical limit. Our numerical experiments confirm our theory and indicate a sharp phase transition for the exact group recovery.
A mass-preserving two-step Lagrange-Galerkin scheme of second order in time for convection-diffusion problems is presented, and convergence with optimal error estimates is proved in the framework of $L^2$-theory. The introduced scheme maintains the advantages of the Lagrange-Galerkin method, i.e., CFL-free robustness for convection-dominated problems and a symmetric and positive coefficient matrix resulting from the discretization. In addition, the scheme conserves the mass on the discrete level if the involved integrals are computed exactly. Unconditional stability and error estimates of second order in time are proved by employing two new key lemmas on the truncation error of the material derivative in conservative form and on a discrete Gronwall inequality for multistep methods. The mass-preserving property is achieved by the Jacobian multiplication technique introduced by Rui and Tabata in 2010, and the accuracy of second order in time is obtained based on the idea of the multistep Galerkin method along characteristics originally introduced by Ewing and Russel in 1981. For the first time step, the mass-preserving scheme of first order in time by Rui and Tabata in 2010 is employed, which is efficient and does not cause any loss of convergence order in the $\ell^\infty(L^2)$- and $\ell^2(H^1_0)$-norms. For the time increment $\Delta t$, the mesh size $h$ and a conforming finite element space of polynomial degree $k$, the convergence order is of $O(\Delta t^2 + h^k)$ in the $\ell^\infty(L^2)\cap \ell^2(H^1_0)$-norm and of $O(\Delta t^2 + h^{k+1})$ in the $\ell^\infty(L^2)$-norm if the duality argument can be employed. Error estimates of $O(\Delta t^{3/2}+h^k)$ in discrete versions of the $L^\infty(H^1_0)$- and $H^1(L^2)$-norm are additionally proved. Numerical results confirm the theoretical convergence orders in one, two and three dimensions.
Within the framework of $ p $-adaptive flux reconstruction, we aim to construct efficient polynomial multigrid ($p$MG) preconditioners for implicit time integration of the Navier--Stokes equations using Jacobian-free Newton--Krylov (JFNK) methods. We hypothesise that in pseudo transient continuation (PTC), as the residual drops, the frequency of error modes that dictates the convergence rate gets higher and higher. We apply nonlinear $p$MG solvers to stiff steady problems at low Mach number ($\mathrm{Ma}=10^{-3}$) to verify our hypothesis. It is demonstrated that once the residual drops by a few orders of magnitude, improved smoothing on intermediate $ p $-sublevels will not only maintain the stability of $ p $MG at large time steps but also improve the convergence rate. For the unsteady Navier--Stokes equations, we elaborate how to construct nonlinear preconditioners using pseudo transient continuation for the matrix-free generalized minimal residual (GMRES) method used in explicit first stage, singly diagonally implicit Runge--Kutta (ESDIRK) methods, and linearly implicit Rosenbrock--Wanner (ROW) methods. Given that at each time step the initial guess in the nonlinear solver is not distant from the converged solution, we recommend a two-level $p\{p_0\text{-}p_0/2\} $ or even $ p\{p_0\text{-}(p_0-1)\} $ $p$-hierarchy for optimal efficiency with a matrix-based smoother on the coarser level based on our hypothesis. It is demonstrated that insufficient smoothing on intermediate $p$-sublevels will deteriorate the performance of $p$MG preconditioner greatly. (See full abstract in the paper.)
This paper is concerned with a numerical solution to the scattering of a time-harmonic electromagnetic wave by a bounded and impenetrable obstacle in three dimensions. The electromagnetic wave propagation is modeled by a boundary value problem of Maxwell's equations in the exterior domain of the obstacle. Based on the Dirichlet-to-Neumann (DtN) operator, which is defined by an infinite series, an exact transparent boundary condition is introduced and the scattering problem is reduced equivalently into a bounded domain. An a posteriori error estimate based adaptive finite element DtN method is developed to solve the discrete variational problem, where the DtN operator is truncated into a sum of finitely many terms. The a posteriori error estimate takes into account both the finite element approximation error and the truncation error of the DtN operator. The latter is shown to decay exponentially with respect to the truncation parameter. Numerical experiments are presented to illustrate the effectiveness of the proposed method.