A C++ library for sensitivity analysis of optimisation problems involving ordinary differential equations (ODEs) enabled by automatic differentiation (AD) and SIMD (Single Instruction, Multiple data) vectorization is presented. The discrete adjoint sensitivity analysis method is implemented for adaptive explicit Runge-Kutta (ERK) methods. Automatic adjoint differentiation (AAD) is employed for efficient evaluations of products of vectors and the Jacobian matrix of the right hand side of the ODE system. This approach avoids the low-level drawbacks of the black box approach of employing AAD on the entire ODE solver and opens the possibility to leverage parallelization. SIMD vectorization is employed to compute the vector-Jacobian products concurrently. We study the performance of other methods and implementations of sensitivity analysis and we find that our algorithm presents a small advantage compared to equivalent existing software.
This paper introduces filtered finite difference methods for numerically solving a dispersive evolution equation with solutions that are highly oscillatory in both space and time. We consider a semiclassically scaled nonlinear Schr\"odinger equation with highly oscillatory initial data in the form of a modulated plane wave. The proposed methods do not need to resolve high-frequency oscillations in both space and time by prohibitively fine grids as would be required by standard finite difference methods. The approach taken here modifies traditional finite difference methods by incorporating appropriate filters. Specifically, we propose the filtered leapfrog and filtered Crank--Nicolson methods, both of which achieve second-order accuracy with time steps and mesh sizes that are not restricted in magnitude by the small semiclassical parameter. Furthermore, the filtered Crank--Nicolson method conserves both the discrete mass and a discrete energy. Numerical experiments illustrate the theoretical results.
The paper considers grad-div stabilized equal-order finite elements (FE) methods for the linearized Navier-Stokes equations. A block triangular preconditioner for the resulting system of algebraic equations is proposed which is closely related to the Augmented Lagrangian (AL) preconditioner. A field-of-values analysis of a preconditioned Krylov subspace method shows convergence bounds that are independent of the mesh parameter variation. Numerical studies support the theory and demonstrate the robustness of the approach also with respect to the viscosity parameter variation, as is typical for AL preconditioners when applied to inf-sup stable FE pairs. The numerical experiments also address the accuracy of grad-div stabilized equal-order FE method for the steady state Navier-Stokes equations.
In this paper we develop numerical analysis for finite element discretization of semilinear elliptic equations with potentially non-Lipschitz nonlinearites. The nonlinearity is essecially assumed to be continuous and monotonically decreasing including the case of non-Lipschitz or even non-H\"older continuous nonlinearities. For a direct discretization (without any regularization) with linear finite elements we derive error estimates with respect to various norms and illustrate these results with a numerical example.
We discuss nonparametric estimation of the trend coefficient in models governed by a stochastic differential equation driven by a multiplicative stochastic volatility.
This work proposes and analyzes an efficient numerical method for solving the nonlinear Schr\"odinger equation with quasiperiodic potential, where the projection method is applied in space to account for the quasiperiodic structure and the Strang splitting method is used in time.While the transfer between spaces of low-dimensional quasiperiodic and high-dimensional periodic functions and its coupling with the nonlinearity of the operator splitting scheme make the analysis more challenging. Meanwhile, compared to conventional numerical analysis of periodic Schr\"odinger systems, many of the tools and theories are not applicable to the quasiperiodic case. We address these issues to prove the spectral accuracy in space and the second-order accuracy in time. Numerical experiments are performed to substantiate the theoretical findings.
Robust and stable high order numerical methods for solving partial differential equations are attractive because they are efficient on modern and next generation hardware architectures. However, the design of provably stable numerical methods for nonlinear hyperbolic conservation laws pose a significant challenge. We present the dual-pairing (DP) and upwind summation-by-parts (SBP) finite difference (FD) framework for accurate and robust numerical approximations of nonlinear conservation laws. The framework has an inbuilt "limiter" whose goal is to detect and effectively resolve regions where the solution is poorly resolved and/or discontinuities are found. The DP SBP FD operators are a dual-pair of backward and forward FD stencils, which together preserve the SBP property. In addition, the DP SBP FD operators are designed to be upwind, that is they come with some innate dissipation everywhere, as opposed to traditional SBP and collocated discontinuous Galerkin spectral element methods which can only induce dissipation through numerical fluxes acting at element interfaces. We combine the DP SBP operators together with skew-symmetric and upwind flux splitting of nonlinear hyperbolic conservation laws. Our semi-discrete approximation is provably entropy-stable for arbitrary nonlinear hyperbolic conservation laws. The framework is high order accurate, provably entropy-stable, convergent, and avoids several pitfalls of current state-of-the-art high order methods. We give specific examples using the in-viscid Burger's equation, nonlinear shallow water equations and compressible Euler equations of gas dynamics. Numerical experiments are presented to verify accuracy and demonstrate the robustness of our numerical framework.
We study the numerical approximation of advection-diffusion equations with highly oscillatory coefficients and possibly dominant advection terms by means of the Multiscale Finite Element Method. The latter method is a now classical, finite element type method that performs a Galerkin approximation on a problem-dependent basis set, itself pre-computed in an offline stage. The approach is implemented here using basis functions that locally resolve both the diffusion and the advection terms. Variants with additional bubble functions and possibly weak inter-element continuity are proposed. Some theoretical arguments and a comprehensive set of numerical experiments allow to investigate and compare the stability and the accuracy of the approaches. The best approach constructed is shown to be adequate for both the diffusion- and advection-dominated regimes, and does not rely on an auxiliary stabilization parameter that would have to be properly adjusted.
We propose a new neural network based large eddy simulation framework for the incompressible Navier-Stokes equations based on the paradigm "discretize first, filter and close next". This leads to full model-data consistency and allows for employing neural closure models in the same environment as where they have been trained. Since the LES discretization error is included in the learning process, the closure models can learn to account for the discretization. Furthermore, we employ a divergence-consistent discrete filter defined through face-averaging and provide novel theoretical and numerical filter analysis. This filter preserves the discrete divergence-free constraint by construction, unlike general discrete filters such as volume-averaging filters. We show that using a divergence-consistent LES formulation coupled with a convolutional neural closure model produces stable and accurate results for both a-priori and a-posteriori training, while a general (divergence-inconsistent) LES model requires a-posteriori training or other stability-enforcing measures.
This paper focuses on the numerical solution of elliptic partial differential equations (PDEs) with Dirichlet and mixed boundary conditions, specifically addressing the challenges arising from irregular domains. Both finite element method (FEM) and finite difference method (FDM), face difficulties in dealing with arbitrary domains. The paper introduces a novel nodal symmetric ghost {method based on a variational formulation}, which combines the advantages of FEM and FDM. The method employs bilinear finite elements on a structured mesh and provides a detailed implementation description. A rigorous a priori convergence rate analysis is also presented. The convergence rates are validated with many numerical experiments, in both one and two space dimensions.
We present a new $hp$-version space-time discontinuous Galerkin (dG) finite element method for the numerical approximation of parabolic evolution equations on general spatial meshes consisting of polygonal/polyhedral (polytopic) elements, giving rise to prismatic space-time elements. A key feature of the proposed method is the use of space-time elemental polynomial bases of \emph{total} degree, say $p$, defined in the physical coordinate system, as opposed to standard dG-time-stepping methods whereby spatial elemental bases are tensorized with temporal basis functions. This approach leads to a fully discrete $hp$-dG scheme using less degrees of freedom for each time step, compared to standard dG time-stepping schemes employing tensorized space-time, with acceptable deterioration of the approximation properties. A second key feature of the new space-time dG method is the incorporation of very general spatial meshes consisting of possibly polygonal/polyhedral elements with \emph{arbitrary} number of faces. A priori error bounds are shown for the proposed method in various norms. An extensive comparison among the new space-time dG method, the (standard) tensorized space-time dG methods, the classical dG-time-stepping, and conforming finite element method in space, is presented in a series of numerical experiments.