A new algorithm is proposed for simultaneously finding multiple solutions of an underdetermined inverse problem. The algorithm was developed for an ODE parameter identification problem in pharmacokinetics for which multiple solutions are of interest. The algorithm proceeds by computing a cluster of solutions simultaneously, and is more efficient than algorithms that compute multiple solutions one-by-one because it fits the Jacobian in a collective way using a least squares approach. It is demonstrated numerically that the algorithm finds accurate solutions that are suitably distributed, guided by a priori information on which part of the solution set is of interest, and that it does so much more efficiently than a baseline Levenberg-Marquardt method that computes solutions one-by-one. It is also demonstrated that the algorithm benefits from improved robustness due to an inherent smoothing provided by the least-squares fitting.
We propose a novel algorithm for the approximation of surface quasi-geostrophic (SQG) flows modeled by a nonlinear partial differential equation coupling transport and fractional diffusion phenomena. The time discretization consists of an explicit strong-stability-preserving three-stage Runge-Kutta method while a flux-corrected-transport (FCT) method coupled with Dunford-Taylor representations of fractional operators is advocated for the space discretization. Standard continuous piecewise linear finite elements are employed, and the algorithm does not have restrictions on the mesh structure or on the computational domain. In the inviscid case, we show that the resulting scheme satisfies a discrete maximum principle property under a standard CFL condition and observe, in practice, its second order accuracy in space. The algorithm successfully approximates several benchmarks with sharp transitions and fine structures typical of SQG flows. In addition, theoretical Kolmogorov energy decay rates are observed on a freely decaying atmospheric turbulence simulation.
Recently, there has been interest in high precision approximations of the first eigenvalue of the Laplace-Beltrami operator on spherical triangles for combinatorial purposes. We compute improved and certified enclosures to these eigenvalues. This is achieved by applying the method of particular solutions in high precision, the enclosure being obtained by a combination of interval arithmetic and Taylor models. The index of the eigenvalue is certified by exploiting the monotonicity of the eigenvalue with respect to the domain. The classically troublesome case of singular corners is handled by combining expansions at all corners and an expansion from an interior point. In particular, this allows us to compute 100 digits of the fundamental eigenvalue for the three-dimensional Kreweras model that has been the object of previous efforts.
In this paper, we develop a family of high order cut discontinuous Galerkin (DG) methods for hyperbolic conservation laws in one space dimension. Ghost penalty stabilization is used to stabilize the scheme for small cut elements. Analysis shows that our proposed methods have similar stability and accuracy properties as the standard DG methods on a regular mesh. We also prove that the cut DG method with piecewise constants in space is total variation diminishing. We use the strong stability preserving Runge-Kutta method for time discretization and the time step is independent on the size of cut element. Numerical examples demonstrate the cut DG methods are high order accurate for smooth problems and perform well for discontinuous problems.
A new second-order method for approximating the compressible Euler equations is introduced. The method preserves all the known invariant domains of the Euler system: positivity of the density, positivity of the internal energy, and the local minimum principle on the specific entropy. The technique combines a first-order, invariant domain preserving, guaranteed maximum speed method using a graph viscosity (GMS-GV1) with an invariant domain violating, but entropy consistent, high-order method. Invariant domain preserving auxiliary states, naturally produced by the GMS-GV1 method, are used to define local bounds for the high-order method, which is then made invariant domain preserving via a convex limiting process. Numerical tests confirm the second-order accuracy of the new GMS-GV2 method in the maximum norm, where the 2 stands for second-order. The proposed convex limiting is generic and can be applied to other approximation techniques and other hyperbolic systems.
In this paper, a new localized radial basis function (RBF) method based on partition of unity (PU) is proposed for solving boundary and initial-boundary value problems. The new method benefits from a direct discretization approach and is called the “direct RBF partition of unity (D-RBF-PU)” method. Thanks to avoiding all derivatives of PU weight functions as well as all lower derivatives of local approximants, the new method is faster and simpler than the standard RBF-PU method. Besides, the discontinuous PU weight functions can now be utilized to develop the method in a more efficient and less expensive way. Alternatively, the new method is an RBF-generated finite difference (RBF-FD) method in a PU setting which is much faster and in some situations more accurate than the original RBF-FD. The polyharmonic splines are used for local approximations, and the error and stability issues are considered. Some numerical experiments on irregular two- and three-dimensional domains, as well as cost comparison tests, are performed to support the theoretical analysis and to show the efficiency of the new method.
We present an implementation of a stage-parallel preconditioner for Radau IIA type fully implicit Runge–Kutta methods, which approximates the inverse of the Runge–Kutta matrix AQ from the Butcher tableau by the lower triangular matrix resulting from an LU decomposition and diagonalizes the system with as many blocks as stages. For the transformed system, we employ a block preconditioner where each block is distributed and solved by a subgroup of processes in parallel. For combination of partial results, we use either a communication pattern resembling Cannon’s algorithm or shared memory. A performance model and a large set of performance studies (including strong-scaling runs with up to 150k processes on 3k compute nodes) conducted for a time-dependent heat problem, using matrix-free finite element methods, indicate that the stage-parallel implementation can reach higher throughputs near the scaling limit. The achievable speedup increases linearly with the number of stages and is bounded by the number of stages. Furthermore, we show that the presented stage-parallel concepts are also applicable to the case that AQ is directly diagonalized, which requires either complex arithmetic or solutions of two-by-two blocks, both exposing about half the parallelism. Alternatively to distributing stages and assigning them to distinct processes, we discuss the possibility of batching operations from different stages together.
An accelerated polynomial expansion scheme to construct the density matrix in quantum mechanical molecular dynamics simulations is proposed. The scheme is based on recursive density matrix expansions, e. g., [A. M. N. Niklasson, Phys. Rev. B, 66 (2002), 155115], which are accelerated by a scale-and-fold technique [E. H. Rubensson, J. Chem. Theory Comput., 7 (2011), pp. 1233-1236]. The acceleration scheme requires interior eigenvalue estimates, which may be expensive and cumbersome to come by. Here we show how such eigenvalue estimates can be extracted from the recursive expansion by a simple and robust procedure at a negligible computational cost. Our method is illustrated with density functional tight-binding Born-Oppenheimer molecular dynamics simulations, where the computational effort is dominated by the density matrix construction. In our analysis we identify two different phases of the recursive polynomial expansion, the conditioning and purification phases, and we show that the acceleration represents an improvement of the conditioning phase, which typically gives a significant reduction of the computational cost.
Localized collocation methods based on radial basis functions (RBFs) for elliptic problems appear to be nonrobust in the presence of Neumann boundary conditions. In this paper, we overcome this issue by formulating the RBF-generated finite difference method in a discrete least squares setting instead. This allows us to prove high-order convergence under node refinement and to numerically verify that the least squares formulation is more accurate and robust than the collocation formulation. The implementation effort for the modified algorithm is comparable to that for the collocation method.