Search

The weighted essentially nonoscillatory (WENO) schemes are popular high-order numerical methods for hyperbolic conservation laws. When dealing with hyperbolic systems, WENO schemes are usually used in cooperation with the local characteristic decomposition, as the componentwise WENO reconstruction/interpolation procedure often produces oscillatory approximations near shocks. In this paper, we investigate local characteristic decomposition-free WENO schemes for a special class of hyperbolic systems endowed with a coordinate system of Riemann invariants. We apply the WENO procedure to the coordinate system of Riemann invariants instead of the local characteristic fields to save the expensive computational cost on local characteristic decomposition but meanwhile maintain the essentially nonoscillatory performance. Due to the nonlinear algebraic relation between the Riemann invariants and conserved variables, it is difficult to obtain the cell averages of Riemann invariants directly from those of conserved variables and vice versa; thus, we do not use the finite volume WENO schemes in this work. The same difficulty is also faced in the traditional Shu-Osher lemma [C.-W. Shu and S. Osher, J. Comput. Phys., 83 (1989), pp. 32-78]-based finite difference schemes, as the computation of fluxes is based on reconstruction as well. Therefore, we adopt the alternative formulation of the finite difference WENO scheme [Y. Jiang, C.-W. Shu, and M. Zhang, SIAM J. Sci. Comput., 35 (2013), pp. A1137-A1160, C.-W. Shu and S. Osher, J. Comput. Phys., 77 (1988), pp. 439-471] in this paper, which is based on interpolation for nodal values. The efficiency and good performance of our method are demonstrated by extensive numerical tests which indicate that the coordinate system of Riemann invariants is a good alternative of local characteristic fields for the WENO procedure. [ABSTRACT FROM AUTHOR]

This paper is focused on the approximation of the Euler equations of compressible fluid dynamics on a staggered mesh. With this aim, the flow parameters are described by the velocity, the density, and the internal energy. The thermodynamic quantities are described on the elements of the mesh, and thus the approximation is only in L2, while the kinematic quantities are globally continuous. The method is general in the sense that the thermodynamic and kinetic parameters are described by an arbitrary degree of polynomials. In practice, the difference between the degrees of the kinematic parameters and the thermodynamic ones is set to 1. The integration in time is done using the forward Euler method but can be extended straightforwardly to higher-order methods. In order to guarantee that the limit solution will be a weak solution of the problem, we introduce a general correction method in the spirit of the Lagrangian staggered method described in [R. Abgrall and S. Tokareva, SIAM J. Sci. Comput., 39 (2017), pp. A2345--A2364; R. Abgrall, K. Lipnikov, N. Morgan, and S. Tokareva, SIAM J. Sci. Comput., 2 (2020), pp. A343--A370; V. A. Dobrev, T. V. Kolev, and R. N. Rieben, SIAM J. Sci. Comput., 34 (2012), pp. B606--B641], and we prove a Lax--Wendroff theorem. The proof is valid for multidimensional versions of the scheme, even though most of the numerical illustrations in this work, on classical benchmark problems, are one-dimensional because we have easy access to the exact solution for comparison. We conclude by explaining that the method is general and can be used in different settings, for example, finite volume or discontinuous Galerkin method, not just the specific one presented in this paper. [ABSTRACT FROM AUTHOR]

Low Reynolds number fluid flows are governed by the Stokes equations. In two dimensions, Stokes flows can be described by two analytic functions, known as Goursat functions. Brubeck and Trefethen [SIAM J. Sci. Comput., 44 (2022), pp. A1205-A1226] recently introduced a lightning Stokes solver that uses rational functions to approximate the Goursat functions in polyg- onal domains. In this paper, we present the LARS algorithm (lightning-AAA rational Stokes) for computing two-dimensional (2D) Stokes flows in domains with smooth boundaries and multiply con- nected domains using lightning and AAA rational approximation [Y. Nakatsukasa, O. Sète, and L. N. Trefethen, SIAM J. Sci. Comput., 40 (2018), pp. A1494-A1522]. After validating our solver against known analytical solutions, we solve a variety of 2D Stokes flow problems with physical and engineering applications. Using these examples, we show rational approximation can now be used to compute 2D Stokes flows in general domains. The computations take less than a second and give solutions with at least 6-digit accuracy. [ABSTRACT FROM AUTHOR]

Entropy inequalities are crucial to the well-posedness of hyperbolic conservation laws, which help to select the physically meaningful one from among the infinite many weak solutions. Recently, several high order discontinuous Galerkin (DG) methods satisfying entropy inequalities were proposed; see [T. Chen and C.-W. Shu, J. Comput. Phys., 345 (2017), pp. 427-461; J. Chan, J. Comput. Phys., 362 (2018), pp. 346-374; T. Chen and C.-W. Shu, CSIAM Trans. Appl. Math., 1 (2020), pp. 1-52] and the references therein. However, high order numerical methods typically generate spurious oscillations in the presence of shock discontinuities. In this paper, we construct a high order entropy stable essentially oscillation-free DG (OFDG) method for hyperbolic conservation laws. With some suitable modification on the high order damping term introduced in [J. Lu, Y. Liu, and C.-W. Shu, SIAM J. Numer. Anal., 59 (2021), pp. 1299-1324; Y. Liu, J. Lu, and C.-W. Shu, SIAM J. Sci. Comput., 44 (2022), pp. A230-A259], we are able to construct an OFDG scheme with dissipative entropy. It is challenging to make the damping term compatible with the current entropy stable DG framework, that is, the damping term should be dissipative for any given entropy function without compromising high order accuracy. The key ingredient is to utilize the convexity of the entropy function and the orthogonality of the projection. Then the proposed method maintains the same properties of conservation, error estimates, and entropy dissipation as the original entropy stable DG method. Extensive numerical experiments are presented to validate the theoretical findings and the capability of controlling spurious oscillations of the proposed algorithm. [ABSTRACT FROM AUTHOR]

The GMRES algorithm of Saad and Schultz [SIAM J. Sci. Stat. Comput., 7 (1986), pp. 856-869] is an iterative method for approximately solving linear systems Ax = b, with initial guess x0 and residual r0 = b Ax0. The algorithm employs the Arnoldi process to generate the Krylov basis vectors (the columns of Vκ). It is well known that this process can be viewed as a QR factorization of the matrix Bκ = [ r0, AVκ] at each iteration. Despite an O (ε)κ (Bκ) loss of orthogonality, for unit roundoff ε and condition number κ, the modified Gram-Schmidt formulation was shown to be backward stable in the seminal paper by Paige et al. [SIAM J. Matrix Anal. Appl., 28 (2006), pp. 264-284]. We present an iterated Gauss-Seidel formulation of the GMRES algorithm (IGS-GMRES) based on the ideas of Ruhe [Linear Algebra Appl., 52 (1983), pp. 591-601] and Świrydowicz et al. [Numer. Linear Algebra Appl., 28 (2020), pp. 1-20]. IGS-GMRES maintains orthogonality to the level O (ε)κ (Bκ) or O (ε), depending on the choice of one or two iterations; for two Gauss-Seidel iterations, the computed Krylov basis vectors remain orthogonal to working accuracy and the smallest singular value of Vκ remains close to one. The resulting GMRES method is thus backward stable. We show that IGS-GMRES can be implemented with only a single synchronization point per iteration, making it relevant to large-scale parallel computing environments. We also demonstrate that, unlike MGS-GMRES, in IGS-GMRES the relative Arnoldi residual corresponding to the computed approximate solution no longer stagnates above machine precision even for highly nonnormal systems. [ABSTRACT FROM AUTHOR]

This paper proposes and analyzes a novel efficient high-order finite volume method for the ideal magnetohydrodynamics (MHD). As a distinctive feature, the method simultaneously preserves two critical physical constraints: a discretely divergence-free (DDF) constraint on the magnetic field and the positivity-preserving (PP) property, which ensures the positivity of density, pressure, and internal energy. To enforce the DDF condition in each cell, we design a new discrete projection approach that projects the reconstructed point values at the cell interface into a DDF space, without using any approximation polynomials. This projection method is highly efficient, easy to implement, and particularly suitable for the high-order finite volume methods that return only the point values (no explicit approximation polynomials) in the reconstruction. Moreover, we also develop a new finite volume framework for constructing provably PP schemes for the ideal MHD system. The framework comprises the discrete projection technique, a suitable approximation to the Godunov--Powell source terms, and a simple PP limiter. We provide rigorous analysis of the PP property of the proposed finite volume method, demonstrating that the DDF condition and the proper approximation to the source terms eliminate the impact of magnetic divergence terms on the PP property. The analysis is challenging due to the internal energy function's nonlinearity and the intricate relationship between the DDF and PP properties. To address these challenges, we adopt the recently developed geometric quasilinearization approach [K. Wu and C.-W. Shu, SIAM Rev., 65 (2023), pp. 1031--1073], which transforms a nonlinear constraint into a family of linear constraints. Finally, we validate the effectiveness of the proposed method through several benchmark and demanding numerical examples. The results demonstrate that the proposed method is robust, accurate, and highly effective, confirming the significance of the proposed DDF projection and PP techniques. [ABSTRACT FROM AUTHOR]

In this paper, we focus our investigation on providing long-term estimates of the Hamiltonian deviation computed along numerical approximations to the solutions of stochastic Hamiltonian systems, of both It\and Stratonovich types. It is well known that the expected Hamiltonian of an It Hamiltonian system with additive noise exhibits a linear drift in time [C. Chen et al., Adv. Comput. Math., 46 (2020), 27], while the Hamiltonian function is conserved along the exact flow of a Stratonovich Hamiltonian system [T. Misawa, Japan J. Indust. Appl. Math., 17 (2000), pp. 119--128; T. Misawa, Math. Probl. Eng., 2010 (2010), 384937]. Here, we focus our attention on providing modified differential equations associated to suitable discretizations for the above problems, by means of weak backward error analysis arguments [K. C. Zygalakis, SIAM J. Sci. Comput., 33 (2011), pp. 102--130]. Then, long-term estimates are provided for both It and Stratonovich Hamiltonian systems, revealing the presence of parasitic terms affecting the overall conservation accuracy. Finally, selected numerical experiments are provided to confirm the theoretical analysis. [ABSTRACT FROM AUTHOR]

This paper builds on the algebraic theory in the companion paper [S. F. McCormick, J. Benzaken, and R. Tamstorf, SIAM J. Sci. Comput., 43 (2021), pp. S392--S419] to obtain discretization-error-accurate solutions for linear elliptic partial differential equations (PDEs) by mixed-precision multigrid solvers. It is often assumed that the achievable accuracy is limited by discretization or algebraic errors. On the contrary, we show that the quantization error incurred by simply storing the matrix in any fixed precision quickly begins to dominate the total error as the discretization is refined. We extend the existing theory to account for these quantization errors and use the resulting bounds to guide the choice of four different precision levels in order to balance quantization, algebraic, and discretization errors in the progressive-precision scheme proposed in the companion paper. A remarkable result is that while iterative refinement is susceptible to quantization errors during the residual and update computations, the V-cycle used to compute the correction in each iteration is much more resilient and continues to work if the system matrices in the hierarchy become indefinite due to quantization. As a result, the V-cycle only requires relatively few bits of precision per level. Based on our findings, we outline a simple way to implement a progressive-precision full multigrid (FMG) solver with minimal overhead and demonstrate as an example that the one-dimensional biharmonic equation can be solved reliably to any desired accuracy using just a few V-cycles when the underlying smoother works well. Additionally, we show that the progressive-precision scheme leads to memory savings of up to 50\% compared to fixed precision. [ABSTRACT FROM AUTHOR]

A general adaptive refinement strategy for solving linear elliptic partial differential equations with random data is proposed and analysed herein. The adaptive strategy extends the a posteriori error estimation framework introduced by Guignard and Nobile [SIAM J. Numer. Anal., 56 (2018), pp. 3121-3143] to cover problems with a nonaffine parametric coefficient dependence. A suboptimal, but nonetheless reliable and convenient implementation of the strategy involves aproximation of the decoupled PDE problems with a common finite element approximation space. Computational results obtained using such a single-level strategy are presented in this paper (part I). Results obtained using a potentially more efficient multilevel approximation strategy, where meshes are individually tailored, will be discussed in part II of this work. The codes used to generate the numerical results are available on GitHub. [ABSTRACT FROM AUTHOR]

This paper proposes some efficient and accurate adaptive two-grid (ATG) finite element algorithms for linear and nonlinear PDEs. The main idea of these algorithms is to utilize the solutions on the kth-level adaptive meshes to find the solutions on the (k + 1)th-level adaptive meshes which are constructed by performing adaptive element bisections on the k th-level adaptive meshes. These algorithms transform nonsymmetric positive definite (non-SPD) PDEs (resp., nonlinear PDEs) into symmetric positive definite (SPD) PDEs (resp., linear PDEs). The proposed algorithms are both accurate and efficient due to the following advantages: They do not need to solve the nonsymmetric or nonlinear systems; the degrees of freedom are very small, they are easily implemented, and the interpolation errors are very small. Next, this paper constructs residual-type a posteriori error estimators which are shown to be reliable and efficient. The key ingredient in proving the efficiency is to establish an upper bound of the oscillation terms, which may not be higher-order terms due to the low regularity of the numerical solution. Furthermore, the convergence of the algorithms is proved when bisection is used for the mesh refinements. Finally, numerical experiments are provided to verify the accuracy and efficiency of the ATG finite element algorithms compared to regular adaptive finite element algorithms and two-grid finite element algorithms [J. Xu, SIAM J. Numer. Anal., 33 (1996), pp. 1759-1777]. [ABSTRACT FROM AUTHOR]

The a posteriori error estimate and greedy algorithms play central roles in the reduced basis method (RBM). In [M. Yano, Comput. Methods Appl. Mech. Engrg., 287 (2015), pp. 290--309; ESAIM Math. Model. Numer. Anal., 50 (2016), pp. 163--185; SIAM J. Sci. Comput., 40 (2018), pp. A388--A420], several versions of RBMs based on exact error certifications and greedy algorithms with spatio-parameter adaptivities are developed. In this paper, with the parametric symmetric coercive elliptic boundary value problem as an example of the primal-dual variational problems satisfying the strong duality, we develop primal-dual RBMs (PD-RBM) with robust true error certifications and discuss three versions of greedy algorithms to balance the finite element error, the exact RB error, and the adaptive mesh refinements. For a class of convex minimization variational problems which has corresponding dual problems satisfying the strong duality, the primal-dual gap between the primal and dual functionals can be used as an a posteriori error estimator. This primal-dual gap error estimator is robust with respect to the parameters of the problem, and it can be used for both mesh refinements of finite element methods and the true RB error certification. With the help of an integration by parts formula, the primal-dual variational theory is developed for the symmetric coercive elliptic boundary value problems with nonhomogeneous boundary conditions by both the conjugate function and Lagrangian theories. A generalized Prager--Synge identity, which is the primal-dual gap error representation for this specific problem, is developed. RBMs for both the primal and dual problems with robust error estimates are developed. The dual variational problem often can be viewed as a constraint optimization problem. In the paper, different from the standard saddlepoint finite element approximation, the dual RBM is treated as a Galerkin projection by constructing RB spaces satisfying the homogeneous constraint. Inspired by the greedy algorithm with spatioparameter adaptivity of [M. Yano, SIAM J. Sci. Comput., 40 (2018), pp. A388--A420], adaptive balanced greedy algorithms with primal-dual finite element and RB error estimators are discussed. Numerical tests are presented to test the PD-RBM with adaptive balanced greedy algorithms. [ABSTRACT FROM AUTHOR]

In this paper, we develop an essentially oscillation-free discontinuous Galerkin method for systems of hyperbolic conservation laws. Based on the standard discontinuous Galerkin method, the numerical damping terms are introduced so as to control the spurious oscillations, similar to the scalar case [J. Lu, Y. Liu, and C.-W. Shu, SIAM J. Numer. Anal., 59 (2021), pp. 1299--1324]. We use both the classical Runge--Kutta method and the modified exponential Runge--Kutta method in time discretization. Particularly, the latter one could avoid additional restrictions of time step size due to the numerical damping. Extensive numerical experiments are shown to demonstrate our algorithm is robust and effective. [ABSTRACT FROM AUTHOR]

In this paper we consider a numerical homogenization technique for curl-curl problems that is based on the framework of the localized orthogonal decomposition and which was proposed in [D. Gallistl, P. Henning, and B. Verfürth, SIAM J. Numer. Anal., 56 (2018), pp. 1570-1596] for problems with essential boundary conditions. The findings of the aforementioned work establish quantitative homogenization results for the time-harmonic Maxwell's equations that hold beyond assumptions of periodicity; however, a practical realization of the approach was left open. In this paper, we transfer the findings from essential boundary conditions to natural boundary conditions, and we demonstrate that the approach yields a computable numerical method. We also investigate how boundary values of the source term can affect the computational complexity and accuracy. Our findings will be supported by various numerical experiments, both in 2D and 3D. [ABSTRACT FROM AUTHOR]

This paper is concerned with direct and inverse scattering by a locally perturbed infinite plane (called a locally rough surface in this paper) on which a Neumann boundary condition is imposed. A novel integral equation formulation is proposed for the direct scattering problem which is defined on a bounded curve (consisting of a bounded part of the infinite plane containing the local perturbation and the lower part of a circle) with two corners and some closed smooth artificial curve. It is a nontrivial extension of our previous work on direct and inverse scattering by a locally rough surface from the Dirichlet boundary condition to the Neumann boundary condition [SIAM J. Appl. Math., 73 (2013), pp. 1811{1829]. For the Dirichlet boundary condition, the integral equation obtained is uniquely solvable in the space of bounded continuous functions on the bounded curve, and it can be solved eficiently by using the Nystrom method with a graded mesh. However, the Neumann condition case leads to an integral equation which is solvable in the space of squarely integrable functions on the bounded curve rather than in the space of bounded continuous functions, making the integral equation very dificult to solve numerically with the classic and eficient Nystrom method. In this paper, we make use of the recursively compressed inverse preconditioning method developed by Helsing to solve the integral equation which is eficient and capable of dealing with large wave numbers. For the inverse problem, it is proved that the locally rough surface is uniquely determined from a knowledge of the far-field pattern corresponding to incident plane waves. Moreover, based on the novel integral equation formulation, a Newton iteration method is developed to reconstruct the locally rough surface from a knowledge of multiple frequency far-field data. Numerical examples are also provided to illustrate that the reconstruction algorithm is stable and accurate even for the case of multiple-scale profiles. [ABSTRACT FROM AUTHOR]

In [McDonald, Pestana, and Wathen, SIAM J. Sci. Comput., 40 (2018), pp. A1012--A1033], a block circulant preconditioner is proposed for all-at-once linear systems arising from evolutionary partial differential equations, in which the preconditioned matrix is proven to be diagonalizable and to have identity-plus-low-rank decomposition in the case of the heat equation. In this paper, we generalize the block circulant preconditioner by introducing a small parameter ϵ>0 into the top-right block of the block circulant preconditioner. The implementation of the generalized preconditioner requires the same computational complexity as that of the block circulant one. Theoretically, we prove that (i) the generalization preserves the diagonalizability and the identity-plus-low-rank decomposition; (ii) all eigenvalues of the new preconditioned matrix are clustered at 1 for sufficiently small\epsilon; (iii) GMRES method for the preconditioned system has a linear convergence rate independent of size of the linear system when\epsilon is taken to be smaller than or comparable to square root of time-step size. Numerical results are reported to confirm the efficiency of the proposed preconditioner and to show that the generalization improves the performance of block circulant preconditioner. [ABSTRACT FROM AUTHOR]

We consider a linear kinetic transport equation under a diffusive scaling that converges to a diffusion equation as the Knudsen number ε → 0. In [S. Boscarino, L. Pareschi, and G. Russo, SIAM J. Sci. Comput., 35 (2013), pp. A22--A51; Z. Peng et al., J. Comput. Phys., 415 (2020), 109485], to achieve the asymptotic preserving (AP) property and unconditional stability in the diffusive regime with ε ≪ 1, numerical schemes are developed based on an additional reformulation of the even-odd or micro-macro decomposed version of the equation. The key of the reformulation is to add a weighted diffusive term on both sides of one equation in the decomposed system. The choice of the weight function, however, is problem-dependent and ad-hoc, and it can affect the performance of numerical simulations. To avoid issues related to the choice of the weight function and still obtain the AP property and unconditional stability in the diffusive regime, we propose in this paper a new family of AP schemes, termed as IMEX-DG-S schemes, directly solving the micro-macro decomposed system without any further reformulation. The main ingredients of the IMEX-DG-S schemes include globally stiffly accurate implicit-explicit (IMEX) Runge--Kutta temporal discretizations with a new IMEX strategy, discontinuous Galerkin spatial discretizations, discrete ordinate methods for the velocity space, and the application of the Schur complement to the algebraic form of the schemes to control the overall computational cost. The AP property of the schemes is shown formally. With an energy type stability analysis applied to the first order scheme, and Fourier type stability analysis applied to the first to third order schemes, we confirm the uniform stability of the methods with respect to\varepsilon and the unconditional stability in the diffusive regime. A series of numerical examples are presented to demonstrate the performance of the new schemes. [ABSTRACT FROM AUTHOR]

The normalized gradient ow, i.e., the gradient ow with discrete normalization (GFDN) introduced in [W. Bao and Q. Du, SIAM J. Sci. Comput., 25 (2004), pp. 1674{1697] or the imaginary time evolution method is one of the most popular techniques for computing the ground states of Bose{Einstein condensates (BECs). In this paper, we revisit the time discretizations for the GFDN and its generalization to the multicomponent BECs. Several widely used time discretizations are demonstrated to be not accurate for computing the ground state solution in the general case, especially for the multicomponent BECs with two or more constraints even for the most accepted linearized backward Euler schemes. More precisely, these schemes usually converge to a solution with an error depending on the time step size. To accurately and efficiently compute the ground state solution of BECs, we propose the gradient ow with Lagrange multiplier (GFLM) method which can be viewed as the modified GFDN by introducing the explicit Lagrange multiplier terms or an approximation of the continuous normalized gradient ow. Through analysis and numerical computation, we clarify that, in order to accurately compute the ground state solution, the GFDN method must be discretized in very special ways in time such as the linearized backward Euler scheme, while the GFLM can be discretized by various schemes and works for multicomponent BEC with multiple constraints. We propose a simple backward-forward Euler scheme to discretize the GFLM, where only a linear system with constant coefficients needs to be solved at each time step. Taking the spin-1 BEC as an example, we demonstrate how the traditional GFDN approach fails in computing the ground state solution of general multicomponent BECs accurately. The proposed GFLM with backward-forward Euler discretization is shown to be accurate and efficient. Extensive numerical results are provided to verify our claims. [ABSTRACT FROM AUTHOR]

In this paper, we present a new regularization scheme for the linearized Poisson{ Boltzmann equation (PBE) which models the electrostatic potential of biomolecules in a solvent. This scheme is based on the splitting of the target potential into the short- and long-range components localized in the molecular region by using the range-separated (RS) tensor format [P. Benner, V. Khoromskaia, and B. N. Khoromskij, SIAM J. Comput., 2 (2018), pp. A1034{A1062] for representation of the discretized multiparticle Dirac delta [B. N. Khoromskij, J. Comput. Phys., 401 (2020), 108998] constituting the highly singular right-hand side in the PBE. From the computational point of view our regularization approach requires only the modification of the right-hand side in the PBE so that it can be implemented within any open-source grid-based software package for solving PBE that already includes some FEM/FDM disretization scheme for elliptic PDE and solver for the arising linear system of equations. The main computational benefits are twofold. First, one applies the chosen PBE solver only for the smooth long-range (regularized) part of the collective potential with the regular right-hand side represented by a low-rank RS tensor with a controllable precision. Thus, we eliminate the numerical treatment of the singularities in the right-hand side and do not change the interface and boundary conditions. And second, the elliptic PDE need not be solved for the singular part in the right-hand side at all, since the short-range part of the target potential of the biomolecule is precomputed independently on a computational grid by simple one-dimensional tensor operations. The total potential is then obtained by adding the numerical solution of the PBE for the smooth long-range part to the directly precomputed tensor representation for the short-range contribution. Numerical tests illustrate that the new regularization scheme, implemented by a simple modification of the right-hand side in the chosen PBE solver, improves the accuracy of the approximate solution on rather coarse grids. The scheme also demonstrates good convergence behavior on a sequence of refined grids. [ABSTRACT FROM AUTHOR]

In this paper, we develop an adaptive multiresolution discontinuous Galerkin (DG) scheme for scalar hyperbolic conservation laws in multidimensions. Compared with previous work for linear hyperbolic equations [W. Guo and Y. Cheng, SIAM J. Sci. Comput., 38 (2016), pp. A3381-A3409, W. Guo and Y. Cheng, SIAM J. Sci. Comput., 39 (2017), pp. A2962-A2992], a class of interpolatory multiwavelets are applied to efficiently compute the nonlinear integrals over elements and edges in DG schemes. The resulting algorithm, therefore, can achieve similar computational complexity as the sparse grid DG method for smooth solutions. Theoretical and numerical studies are performed taking into consideration the accuracy and stability with regard to the choice of the interpolatory multiwavelets. Artificial viscosity is added to capture the shock and only acts on the leaf elements taking advantage of the multiresolution representation. Adaptivity is realized by auto error thresholding based on hierarchical surplus. Accuracy and robustness are demonstrated by several numerical tests. [ABSTRACT FROM AUTHOR]

The Strang splitting method, formally of order two, can suffer from order reduction when applied to semilinear parabolic problems with inhomogeneous boundary conditions. The recent work [L. Einkemmer and A. Ostermann, SIAM J. Sci. Comput., 37, 2015; SIAM J. Sci. Comput., 38, 2016] introduces a modification of the method to avoid the reduction of order based on the nonlinearity. In this paper we introduce a new correction constructed directly from the flow of the nonlinearity and which requires no evaluation of the source term or its derivatives. The goal is twofold. One, this new modification requires only one evaluation of the diffusion flow and one evaluation of the source term flow at each step of the algorithm and it reduces the computational effort to construct the correction. Second, numerical experiments suggest it is well suited in the case where the nonlinearity is stiff. We provide a convergence analysis of the method for a smooth nonlinearity and perform numerical experiments to illustrate the performances of the new approach. [ABSTRACT FROM AUTHOR]

We present extensions to hierarchical probing, a method developed in [A. Stathopoulos, J. Laeuchli, and K. Orginos, SIAM J. Sci. Comput., 35 (2013), pp. S299-S322] to reduce the variance of the Monte Carlo estimation of the trace or the diagonal of the inverse of a large, sparse matrix. In that context, probing is a method to determine the largest-in-magnitude elements of the matrix inverse and then annihilate their contributions to the variance by solving linear systems with appropriate probing vectors. It typically involves coloring the graph of An, since this matches the sparsity structure of a polynomial approximation to A-1. This is equivalent to distance-n coloring of A, i.e., determining which nodes are connected to one other at distance ≤ n. For matrices that display a Green's function decay, n is small, which reduces the number of linear systems to be solved. Our hierarchical probing method was developed for matrices with a lattice structure, where distancen coloring and the generation of probing vectors can be performed far more efficiently and in a way so that earlier vectors are subsets of vectors generated later in the process, meaning that it is simple to continue probing if additional accuracy is needed. However, this method worked only on lattices with dimension lengths that were powers of two. In this paper we extend the method to work on lattices of arbitrary dimension lengths, which is theoretically more challenging. Additionally, we expand the idea to a multilevel, hierarchical probing heuristic for matrices with any undirected graph structure that matches the performance of classical probing but with tractable memory requirements. [ABSTRACT FROM AUTHOR]

Optimization-based samplers such as randomize-then-optimize (RTO) [J. M. Bardsley et al., SIAM J. Sci. Comput., 36 (2014), pp. A1895--A1910] provide an efficient and parallellizable approach to solving large-scale Bayesian inverse problems. These methods solve randomly perturbed optimization problems to draw samples from an approximate posterior distribution. ''Correcting"" these samples, either by Metropolization or importance sampling, enables characterization of the original posterior distribution. This paper focuses on the scalability of RTO to problems with highor infinite-dimensional parameters. In particular, we introduce a new subspace strategy to reformulate RTO. For problems with intrinsic low-rank structures, this subspace acceleration makes the computational complexity of RTO scale linearly with the parameter dimension. Furthermore, this subspace perspective suggests a natural extension of RTO to a function space setting. We thus formalize a function space version of RTO and establish sufficient conditions for it to produce a valid Metropolis--Hastings proposal, yielding dimension-independent sampling performance. Numerical examples corroborate the dimension independence of RTO and demonstrate sampling performance that is also robust to small observational noise. [ABSTRACT FROM AUTHOR]

Several reduced order models have been successfully developed for nonlinear dynamical systems. To achieve a considerable speed-up, a hyper-reduction step is needed to reduce the computational complexity due to nonlinear terms. Many hyper-reduction techniques require the construction of nonlinear term basis, which introduces a computationally expensive offline phase. A novel way of constructing nonlinear term basis within the hyper-reduction process is introduced. In contrast to the traditional hyper-reduction techniques where the collection of nonlinear term snapshots is required, the SNS method avoids collecting the nonlinear term snapshots. Instead, it uses the solution snapshots that are used for building a solution basis, which enables avoiding an extra data compression of nonlinear term snapshots. As a result, the SNS method provides a more efficient offline strategy than the traditional model order reduction techniques, such as the DEIM, GNAT, and ST-GNAT methods. The SNS method is theoretically justified by the conforming subspace condition and the subspace inclusion relation. It is useful for model order reduction of large-scale nonlinear dynamical problems to reduce the offline cost. It is especially useful for ST-GNAT that has shown promising results, such as a good accuracy with a considerable online speed-up for hyperbolic problems in a recent paper by Choi and Carlberg [SIAM J. Sci. Comput., 41 (2019), pp. A26--A58], because ST-GNAT involves an expensive offline cost related to collecting nonlinear term snapshots. Error analysis for the SNS method is presented. Numerical results support that the accuracy of the solution from the SNS method is comparable to the traditional methods and a considerable speed-up (i.e., a factor of two to a hundred) is achieved in the offline phase. [ABSTRACT FROM AUTHOR]

It is shown in [SIAM J. Sci. Comput., 39 (2017), pp. B424-B441] that free-form curves used in computer-aided geometric design can usually be represented as the solutions of linear differential systems and that points and derivatives on the curves can be evaluated dynamically by solving the differential systems numerically. In this paper we present an even more robust and efficient algorithm for dynamic evaluation of exponential polynomial curves and surfaces. Based on properties that spaces spanned by general exponential polynomials are translation invariant and polynomial spaces are invariant with respect to a linear transformation of the parameter, the transformation matrices between bases with or without translated or linearly transformed parameters are explicitly computed. Points on curves or surfaces with equal or changing parameter steps can then be evaluated dynamically from a start point using a precomputed matrix. Like former dynamic evaluation algorithms, the newly proposed approach needs only arithmetic operations for evaluating exponential polynomial curves and surfaces. Unlike conventional numerical methods that solve a linear differential system, the new method can give robust and accurate evaluation results for any chosen parameter steps. The basis transformation technique also enables dynamic evaluation of polynomial curves with changing parameter steps using a constant matrix, which reduces time costs significantly than computing each point individually by classical algorithms. [ABSTRACT FROM AUTHOR]

In [L. Li, M. Ni, and W. Zheng, SIAM J. Sci. Comput., 41 (2019), pp. B796-B815] a charge-conservative finite element method is proposed for solving inductionless and incompressible magnetohydrodynamic (MHD) equations. The purpose of this paper is to propose a robust solver for the discrete problem. Using the framework of field-of-values-equivalence, we first study the preconditioned Krylov space method for the continuous problem in the setting of Hilbert spaces. The algebraic preconditioner for the discrete problem is then obtained by representing the preconditioner for the continuous problem in finite element spaces. By three numerical examples, the optimality of the solver to the number of unknowns is demonstrated for both stationary and time-dependent MHD problems. [ABSTRACT FROM AUTHOR]