GPU Accelerated Discontinuous Galerkin Methods for Shallow Water Equations

We discuss the development, verification, and performance of a GPU accelerated discontinuous Galerkin method for the solutions of two dimensional nonlinear shallow water equations. The shallow water equations are hyperbolic partial differential equations and are widely used in the simulation of tsunami wave propagations. Our algorithms are tailored to take advantage of the single instruction multiple data (SIMD) architecture of graphic processing units. The time integration is accelerated by local time stepping based on a multi-rate Adams-Bashforth scheme. A total variational bounded limiter is adopted for nonlinear stability of the numerical scheme. This limiter is coupled with a mass and momentum conserving positivity preserving limiter for the special treatment of a dry or partially wet element in the triangulation. Accuracy, robustness and performance are demonstrated with the aid of test cases. We compare the performance of the kernels expressed in a portable threading language OCCA, when cross compiled with OpenCL, CUDA, and OpenMP at runtime.Comment: 26 pages, 51 figure

Numerical methods and applications for reduced models of blood flow

The human cardiovascular system is a vastly complex collection of interacting components, including vessels, organ systems, valves, regulatory mechanisms, microcirculations, remodeling tissue, and electrophysiological signals. Experimental, mathematical, and computational research efforts have explored various hemodynamic questions; the scope of this literature is a testament to the intricate nature of cardiovascular physiology. In this work, we focus on computational modeling of blood flow in the major vessels of the human body. We consider theoretical questions related to the numerical approximation of reduced models for blood flow, posed as nonlinear hyperbolic systems in one space dimension. Further, we apply this modeling framework to abnormal physiologies resulting from surgical intervention in patients with congenital heart defects. This thesis contains three main parts: (i) a discussion of the implementation and analysis for numerical discretizations of reduced models for blood flow, (ii) an investigation of solutions to different classes of models in the realm of smooth and discontinuous solutions, and (iii) an application of these models within a multiscale framework for simulating flow in patients with hypoplastic left heart syndrome. The two numerical discretizations studied in this thesis are a characteristics-based method for approximating the Riemann-invariants of reduced blood flow models, and a discontinuous Galerkin scheme for approximating solutions to the reduced models directly. A priori error estimates are derived in particular cases for both methods. Further, two classes of hyperbolic systems for blood flow, namely the mass-momentum and the mass-velocity formulations, are systematically compared with each numerical method and physiologically relevant networks of vessels and boundary conditions. Lastly, closed loop vessel network models of various Fontan physiologies are constructed. Arterial and venous trees are built from networks of one-dimensional vessels while the heart, valves, vessel junctions, and organ beds are modeled by systems of algebraic and ordinary differential equations

Implicit-explicit multistep formulations for finite element discretisations using continuous interior penalty

We consider a finite element method with symmetric stabilisation for the discretisation of the transient convection–diffusion equation. For the time-discretisation we consider either the second order backwards differentiation formula or the Crank–Nicolson method. Both the convection term and the associated stabilisation are treated explicitly using an extrapolated approximate solution. We prove stability of the method and the t2+hp+12 error estimates for the L2-norm under either the standard hyperbolic CFL condition, when piecewise affine (p=1) approximation is used, or in the case of finite element approximation of order p≥1, a stronger, so-called 4/3-CFL, i.e. t≤Ch4/3. The theory is illustrated with some numerical examples

The Discontinuous Galerkin Finite Element Method for Ordinary Differential Equations

We present an analysis of the discontinuous Galerkin (DG) finite element method for nonlinear ordinary differential equations (ODEs). We prove that the DG solution is $(p + 1)$th order convergent in the $L^2$-norm, when the space of piecewise polynomials of degree $p$ is used. A $(2p+1)$th order superconvergence rate of the DG approximation at the downwind point of each element is obtained under quasi-uniform meshes. Moreover, we prove that the DG solution is superconvergent with order $p+2$ to a particular projection of the exact solution. The superconvergence results are used to show that the leading term of the DG error is proportional to the $(p + 1)$-degree right Radau polynomial. These results allow us to develop a residual-based a posteriori error estimator which is computationally simple, efficient, and asymptotically exact. The proposed a posteriori error estimator is proved to converge to the actual error in the $L^2$-norm with order $p+2$. Computational results indicate that the theoretical orders of convergence are optimal. Finally, a local adaptive mesh refinement procedure that makes use of our local a posteriori error estimate is also presented. Several numerical examples are provided to illustrate the global superconvergence results and the convergence of the proposed estimator under mesh refinement

Analysis of \u3ci\u3ea posteriori\u3c/i\u3e error estimates of the discontinuous Galerkin method for nonlinear ordinary differential equations

We develop and analyze a new residual-based a posteriori error estimator for the discontinuous Galerkin (DG) method for nonlinear ordinary differential equations (ODEs). The a posteriori DG error estimator under investigation is computationally simple, efficient, and asymptotically exact. It is obtained by solving a local residual problem with no boundary condition on each element. We first prove that the DG solution exhibits an optimal O(hp+1) convergence rate in the L2-norm when p-degree piece-wise polynomials with p ≥1 are used. We further prove that the DG solution is O(h2p+1) superconvergent at the downwind points. We use these results to prove that the p-degree DG solution is O(hp+2) super close to a particular projection of the exact solution. This superconvergence result allows us to show that the true error can be divided into a significant part and a less significant part. The significant part of the discretization error for the DG solution is proportional to the (p +1)-degree right Radau polynomial and the less significant part converges at O(hp+2) rate in the L2-norm. Numerical experiments demonstrate that the theoretical rates are optimal. Based on the global superconvergent approximations, we construct asymptotically exact a posteriori error estimates and prove that they converge to the true errors in the L2-norm under mesh refinement. The order of convergence is proved to be p +2. Finally, we prove that the global effectivity index in the L2-norm converges to unity at O(h)rate. Several numerical examples are provided to illustrate the global superconvergence results and the convergence of the proposed estimator under mesh refinement. A local adaptive procedure that makes use of our local a posteriori error estimate is also presented

Nonstandard Finite Element Methods

[no abstract available

Uncertainty Quantification by MLMC and Local Time-stepping For Wave Propagation

Because of their robustness, efficiency and non-intrusiveness, Monte Carlo methods are probably the most popular approach in uncertainty quantification to computing expected values of quantities of interest (QoIs). Multilevel Monte Carlo (MLMC) methods significantly reduce the computational cost by distributing the sampling across a hierarchy of discretizations and allocating most samples to the coarser grids. For time dependent problems, spatial coarsening typically entails an increased time-step. Geometric constraints, however, may impede uniform coarsening thereby forcing some elements to remain small across all levels. If explicit time-stepping is used, the time-step will then be dictated by the smallest element on each level for numerical stability. Hence, the increasingly stringent CFL condition on the time-step on coarser levels significantly reduces the advantages of the multilevel approach. By adapting the time-step to the locally refined elements on each level, local time-stepping (LTS) methods permit to restore the efficiency of MLMC methods even in the presence of complex geometry without sacrificing the explicitness and inherent parallelism

Uncertainty quantification by multilevel Monte Carlo and local time-stepping

Because of their robustness, efficiency, and non intrusiveness, Monte Carlo methods are probablythe most popular approach in uncertainty quantification for computing expected values of quantitiesof interest. Multilevel Monte Carlo (MLMC) methods significantly reduce the computational costby distributing the sampling across a hierarchy of discretizations and allocating most samples tothe coarser grids. For time dependent problems, spatial coarsening typically entails an increasedtime step. Geometric constraints, however, may impede uniform coarsening thereby forcing someelements to remain small across all levels. If explicit time-stepping is used, the time step will thenbe dictated by the smallest element on each level for numerical stability. Hence, the increasinglystringent CFL condition on the time step on coarser levels significantly reduces the advantages of themultilevel approach. To overcome that bottleneck we propose to combine the multilevel approach ofMLMC with local time-stepping. By adapting the time step to the locally refined elements on eachlevel, the efficiency of MLMC methods is restored even in the presence of complex geometry withoutsacrificing the explicitness and inherent parallelism. In a careful cost comparison, we quantify thereduction in computational cost for local refinement either inside a small fixed region or towards areentrant corner

Weighted Error Estimates for Transient Transport Problems Discretized Using Continuous Finite Elements with Interior Penalty Stabilization on the Gradient Jumps

In this paper we consider the semi-discretization in space of a first order scalar transport equation. For the space discretization we use standard continuous finite elements with a stabilization consisting of a penalty on the jump of the gradient over element faces. We recall some global error estimates for smooth and rough solutions and then prove a new local error estimate for the transient linear transport equation. In particular we show that for the stabilized method the effect of non-smooth features in the solution decay exponentially from the space time zone where the solution is rough so that smooth features will be transported unperturbed. Locally the L2-norm error converges with the expected order O(hk+12), if the exact solution is locally smooth. We then illustrate the results numerically. In particular we show the good local accuracy in the smooth zone of the stabilized method and that the standard Galerkin fails to approximate a solution that is smooth at the final time if underresolved features have been present in the solution at some time during the evolution

Solution strategies for nonlinear conservation laws

Nonlinear conservation laws form the basis for models for a wide range of physical phenomena. Finding an optimal strategy for solving these problems can be challenging, and a good strategy for one problem may fail spectacularly for others. As different problems have different challenging features, exploiting knowledge about the problem structure is a key factor in achieving an efficient solution strategy. Most strategies found in literature for solving nonlinear problems involve a linearization step, usually using Newton's method, which replaces the original nonlinear problem by an iteration process consisting of a series of linear problems. A large effort is then spent on finding a good strategy for solving these linear problems. This involves choosing suitable preconditioners and linear solvers. This approach is in many cases a good choice and a multitude of different methods have been developed. However, the linearization step to some degree involves a loss of information about the original problem. This is not necessarily critical, but in many cases the structure of the nonlinear problem can be exploited to a larger extent than what is possible when working solely on the linearized problem. This may involve knowledge about dominating physical processes and specifically on whether a process is near equilibrium. By using nonlinear preconditioning techniques developed in recent years, certain attractive features such as automatic localization of computations to parts of the problem domain with the highest degree of nonlinearities arise. In the present work, these methods are further refined to obtain a framework for nonlinear preconditioning that also takes into account equilibrium information. This framework is developed mainly in the context of porous media, but in a general manner, allowing for application to a wide range of problems. A scalability study shows that the method is scalable for challenging two-phase flow problems. It is also demonstrated for nonlinear elasticity problems. Some models arising from nonlinear conservation laws are best solved using completely different strategies than the approach outlined above. One such example can be found in the field of surface gravity waves. For special types of nonlinear waves, such as solitary waves and undular bores, the well-known Korteweg-de Vries (KdV) equation has been shown to be a suitable model. This equation has many interesting properties not typical of nonlinear equations which may be exploited in the solver, and strategies usually reserved to linear problems may be applied. In this work includes a comparative study of two discretization methods with highly different properties for this equation