Multi-Scale Deep Neural Network (MscaleDNN) Methods for Oscillatory Stokes Flows in Complex Domains

In this paper, we study a multi-scale deep neural network (MscaleDNN) as a meshless numerical method for computing oscillatory Stokes flows in complex domains. The MscaleDNN employs a multi-scale structure in the design of its DNN using radial scalings to convert the approximation of high frequency components of the highly oscillatory Stokes solution to one of lower frequencies. The MscaleDNN solution to the Stokes problem is obtained by minimizing a loss function in terms of $L^2$ norm of the residual of the Stokes equation. Three forms of loss functions are investigated based on vorticity-velocity-pressure, velocity-stress-pressure, and velocity-gradient of velocity-pressure formulations of the Stokes equation. We first conduct a systematic study of the MscaleDNN methods with various loss functions on the Kovasznay flow in comparison with normal fully connected DNNs. Then, Stokes flows with highly oscillatory solutions in a 2-D domain with six randomly placed holes are simulated by the MscaleDNN. The results show that MscaleDNN has faster convergence and consistent error decays in the simulation of Kovasznay flow for all three tested loss functions. More importantly, the MscaleDNN is capable of learning highly oscillatory solutions when the normal DNNs fail to converge.     

Adaptive Finite Element Modeling Techniques for the Poisson-Boltzmann Equation

We consider the design of an effective and reliable adaptive finite element method (AFEM) for the nonlinear Poisson-Boltzmann equation (PBE). We first examine the two-term regularization technique for the continuous problem recently proposed by Chen, Holst and Xu based on the removal of the singular electrostatic potential inside biomolecules; this technique made possible the development of the first complete solution and approximation theory for the Poisson-Boltzmann equation, the first provably convergent discretization and also allowed for the development of a provably convergent AFEM. However, in practical implementation, this two-term regularization exhibits numerical instability. Therefore, we examine a variation of this regularization technique which can be shown to be less susceptible to such instability. We establish a priori estimates and other basic results for the continuous regularized problem, as well as for Galerkin finite element approximations. We show that the new approach produces regularized continuous and discrete problems with the same mathematical advantages of the original regularization. We then design an AFEM scheme for the new regularized problem and show that the resulting AFEM scheme is accurate and reliable, by proving a contraction result for the error. This result, which is one of the first results of this type for nonlinear elliptic problems, is based on usingcontinuous and discrete a priori L^∞ estimates. To provide a high-quality geometric model as input to the AFEM algorithm, we also describe a class of feature-preserving adaptive mesh generation algorithms designed specifically for constructing meshes of biomolecular structures, based on the intrinsic local structure tensor of the molecular surface. All of the algorithms described in the article are implemented in the Finite Element Toolkit (FETK), developed and maintained at UCSD. The stability advantages of the new regularization scheme are demonstrated with FETK through comparisons with the original regularization approach for a model problem. The convergence and accuracy of the overall AFEM algorithm is also illustrated by numerical approximation of electrostatic solvation energy for an insulin protein.     

A Kernel-Free Boundary Integral Method for Variable Coefficients Elliptic PDEs

This work proposes a generalized boundary integral method for variable coefficients elliptic partial differential equations (PDEs), including both boundary value and interface problems. The method is kernel-free in the sense that there is no need to know analytical expressions for kernels of the boundary and volume integrals in the solution of boundary integral equations. Evaluation of a boundary or volume integral is replaced with interpolation of a Cartesian grid based solution, which satisfies an equivalent discrete interface problem, while the interface problem is solved by a fast solver in the Cartesian grid. The computational work involved with the generalized boundary integral method is essentially linearly proportional to the number of grid nodes in the domain. This paper gives implementation details for a second-order version of the kernel-free boundary integral method in two space dimensions and presents numerical experiments to demonstrate the efficiency and accuracy of the method for both boundary value and interface problems. The interface problems demonstrated include those with piecewise constant and large-ratio coefficients and the heterogeneous interface problem, where the elliptic PDEs on two sides of the interface are of different types.     

A Fourth-Order Kernel-Free Boundary Integral Method for Interface Problems

This paper presents a fourth-order Cartesian grid based boundary integral method (BIM) for heterogeneous interface problems in two and three dimensional space, where the problem interfaces are irregular and can be explicitly given by parametric curves or implicitly defined by level set functions. The method reformulates the governing equation with interface conditions into boundary integral equations (BIEs) and reinterprets the involved integrals as solutions to some simple interface problems in an extended regular region. Solution of the simple equivalent interface problems for integral evaluation relies on a fourth-order finite difference method with an FFT-based fast elliptic solver. The structure of the coefficient matrix is preserved even with the existence of the interface. In the whole calculation process, analytical expressions of Green’s functions are never determined, formulated or computed. This is the novelty of the proposed kernel-free boundary integral (KFBI) method. Numerical experiments in both two and three dimensions are shown to demonstrate the algorithm efficiency and solution accuracy even for problems with a large diffusion coefficient ratio.     

Convergent Adaptive Finite Element Method Based on Centroidal Voronoi Tessellations and Superconvergence

We present a novel adaptive finite element method (AFEM) for elliptic equations which is based upon the Centroidal Voronoi Tessellation (CVT) and superconvergent gradient recovery. The constructions of CVT and its dual Centroidal Voronoi Delaunay Triangulation (CVDT) are facilitated by a localized Lloyd iteration to produce almost equilateral two dimensional meshes. Working with finite element solutions on such high quality triangulations, superconvergent recovery methods become particularly effective so that asymptotically exact a posteriori error estimations can be obtained. Through a seamless integration of these techniques, a convergent adaptive procedure is developed. As demonstrated by the numerical examples, the new AFEM is capable of solving a variety of model problems and has great potential in practical applications.     

A priori error estimates for optimal control problems with state and control constraints

The Galerkin spectral method is used to approximate elliptic optimal control problems with integral state and control constraints in this paper. We derive the optimality conditions for the control problems by E. Casas' results, then a priori error estimates are provided. Some numerical examples are investigated to confirm the theoretical results and to show the efficiency of the proposed method.     

An Augmented Lagrangian Deep Learning Method for Variational Problems with Essential Boundary Conditions

This paper is concerned with a novel deep learning method for variational problems with essential boundary conditions. To this end, we first reformulate the original problem into a minimax problem corresponding to a feasible augmented Lagrangian, which can be solved by the augmented Lagrangian method in an infinite dimensional setting. Based on this, by expressing the primal and dual variables with two individual deep neural network functions, we present an augmented Lagrangian deep learning method for which the parameters are trained by the stochastic optimization method together with a projection technique. Compared to the traditional penalty method, the new method admits two main advantages: i) the choice of the penalty parameter is flexible and robust, and ii) the numerical solution is more accurate in the same magnitude of computational cost. As typical applications, we apply the new approach to solve elliptic problems and (nonlinear) eigenvalue problems with essential boundary conditions, and numerical experiments are presented to show the effectiveness of the new method.     

Parareal-Richardson Algorithm for Solving Nonlinear ODEs and PDEs

The parareal algorithm, proposed firstly by Lions et al. [J. L. Lions, Y. Maday, and G. Turinici, A "parareal" in time discretization of PDE's, C.R. Acad. Sci. Paris Sér. I Math., 332 (2001), pp. 661-668], is an effective algorithm to solve the time-dependent problems parallel in time. This algorithm has received much interest from many researchers in the past years. We present in this paper a new variant of the parareal algorithm, which is derived by combining the original parareal algorithm and the Richardson extrapolation, for the numerical solution of the nonlinear ODEs and PDEs. Several nonlinear problems are tested to show the advantage of the new algorithm. The accuracy of the obtained numerical solution is compared with that of its original version (i.e., the parareal algorithm based on the same numerical method).     

Operator Splitting for Three-Phase Flow in Heterogeneous Porous Media

We describe an operator splitting technique based on physics rather than on dimension for the numerical solution of a nonlinear system of partial differential equations which models three-phase flow through heterogeneous porous media. The model for three-phase flow considered in this work takes into account capillary forces, general relations for the relative permeability functions and variable porosity and permeability fields. In our numerical procedure a high resolution, nonoscillatory, second order, conservative central difference scheme is used for the approximation of the nonlinear system of hyperbolic conservation laws modeling the convective transport of the fluid phases. This scheme is combined with locally conservative mixed finite elements for the numerical solution of the parabolic and elliptic problems associated with the diffusive transport of fluid phases and the pressure-velocity problem. This numerical procedure has been used to investigate the existence and stability of nonclassical shock waves (called transitional or undercompressive shock waves) in two-dimensional heterogeneous flows, thereby extending previous results for one-dimensional flow problems. Numerical experiments indicate that the operator splitting technique discussed here leads to computational efficiency and accurate numerical results.     

Two-Grid Method for Miscible Displacement Problem by Mixed Finite Element Methods and Mixed Finite Element Method of Characteristics

The miscible displacement of one incompressible fluid by another in a porous medium is governed by a system of two equations. One is elliptic form equation for the pressure and the other is parabolic form equation for the concentration of one of the fluids. Since only the velocity and not the pressure appears explicitly in the concentration equation, we use a mixed finite element method for the approximation of the pressure equation and mixed finite element method with characteristics for the concentration equation. To linearize the mixed-method equations, we use a two-grid algorithm based on the Newton iteration method for this full discrete scheme problems. First, we solve the original nonlinear equations on the coarse grid, then, we solve the linearized problem on the fine grid used Newton iteration once. It is shown that the coarse grid can be much coarser than the fine grid and achieve asymptotically optimal approximation as long as the mesh sizes satisfy $h=H^2$ in this paper. Finally, numerical experiment indicates that two-grid algorithm is very effective.     

A Least-Squares/Fictitious Domain Method for Linear Elliptic Problems with Robin Boundary Conditions

In this article, we discuss a least-squares/fictitious domain method for the solution of linear elliptic boundary value problems with Robin boundary conditions. Let Ω and ω be two bounded domains of R^dsuch that ω⊂Ω. For a linear elliptic problem in Ω\ω with Robin boundary condition on the boundary γ of ω, our goal here is to develop a fictitious domain method where one solves a variant of the original problem on the full Ω, followed by a well-chosen correction over ω. This method is of the virtual control type and relies on a least-squares formulation making the problem solvable by a conjugate gradient algorithm operating in a well chosen control space. Numerical results obtained when applying our method to the solution of two-dimensional elliptic and parabolic problems are given; they suggest optimal order of convergence.     

An Energy Regularization for Cauchy Problems of Laplace Equation in Annulus Domain

Detecting corrosion by electrical field can be modeled by a Cauchy problem of Laplace equation in annulus domain under the assumption that the thickness of the pipe is relatively small compared with the radius of the pipe. The interior surface of the pipe is inaccessible and the nondestructive detection is solely based on measurements from the outer layer. The Cauchy problem for an elliptic equation is a typical ill-posed problem whose solution does not depend continuously on the boundary data. In this work, we assume that the measurements are available on the whole outer boundary on an annulus domain. By imposing reasonable assumptions, the theoretical goal here is to derive the stabilities of the Cauchy solutions and an energy regularization method. Relationship between the proposed energy regularization method and the Tikhonov regularization with Morozov principle is also given. A novel numerical algorithm is proposed and numerical examples are given.     

Enforcing the Discrete Maximum Principle for Linear Finite Element Solutions of Second-Order Elliptic Problems

The maximum principle is a basic qualitative property of the solution of second-order elliptic boundary value problems. The preservation of the qualitative characteristics, such as the maximum principle, in discrete model is one of the key requirements. It is well known that standard linear finite element solution does not satisfy maximum principle on general triangular meshes in 2D. In this paper we consider how to enforce discrete maximum principle for linear finite element solutions for the linear second-order self-adjoint elliptic equation. First approach is based on repair technique, which is a posteriori correction of the discrete solution. Second method is based on constrained optimization. Numerical tests that include anisotropic cases demonstrate how our method works for problems for which the standard finite element methods produce numerical solutions that violate the discrete maximum principle.     

A Parallel Second Order Cartesian Method for Elliptic Interface Problems

We present a parallel Cartesian method to solve elliptic problems with complex immersed interfaces. This method is based on a finite-difference scheme and is second-order accurate in the whole domain. The originality of the method lies in the use of additional unknowns located on the interface, allowing to express straightforwardly the interface transmission conditions. We describe the method and the details of its parallelization performed with the PETSc library. Then we present numerical validations in two dimensions, assorted with comparisons to other related methods, and a numerical study of the parallelized method.     

Development of a Volume of Fluid Method for Computing Interfacial Incompressible Fluid Flows

This study is aimed to develop a volume of fluid (VOF) method to capture the free surface flow. The incompressible two-phase flow is computed by second-order Adams-Bashforth algorithm with a uniform staggered Cartesian grid arrangement. The tangent of hyperbola for interface capturing (THINC) scheme and weighted linear interface calculation (WLIC) based geometrical reconstruction procedure have been implemented in the operator-splitting method for the VOF method. The proposed VOF method preserves mass well, and the interface normal vector can be easily estimated from the level set (LS) function. The LS function, which is a continuous signed distance function around the interface, is represented by solving the re-initialization equation. Numerical results using the present scheme are compared with experimental data and other numerical results in the Rayleigh-Taylor instability, dam-break flow, travelling solitary wave, Kelvin-Helmholtz instability, rising bubble and merging bubble problems. We also present numerical results in detail between computations made with the proposed VOF method and computations made with the conventional LS method.     

dsoa: The implementation of a dynamic system optimization algorithm

This paper describes the ANSI C/C++ computer program dsoa , which implements an algorithm for the approximate solution of dynamics system optimization problems. The algorithm is a direct method that can be applied to the optimization of dynamic systems described by index‐1 differential‐algebraic equations (DAEs). The types of problems considered include optimal control problems and parameter identification problems. The numerical techniques are employed to transform the dynamic system optimization problem into a parameter optimization problem by: (i) parameterizing the control input as piecewise constant on a fixed mesh, and (ii) approximating the DAEs using a linearly implicit Runge‐Kutta method. The resultant nonlinear programming (NLP) problem is solved via a sequential quadratic programming technique. The program dsoa is evaluated using 83 nontrivial optimal control problems that have appeared in the literature. Here we compare the performance of the algorithm using two different NLP problem solvers, and two techniques for computing the derivatives of the functions that define the problem. Copyright © 2009 John Wiley & Sons, Ltd.     

Modeling 3D Magma Dynamics Using a Discontinuous Galerkin Method

Discontinuous Galerkin (DG) and matrix-free finite element methods with a novel projective pressure estimation are combined to enable the numerical modeling of magma dynamics in 2D and 3D using the library deal.II. The physical model is an advection-reaction type system consisting of two hyperbolic equations to evolve porosity and soluble mineral abundance at local chemical equilibrium and one elliptic equation to recover global pressure. A combination of a discontinuous Galerkin method for the advection equations and a finite element method for the elliptic equation provide a robust and efficient solution to the channel regime problems of the physical system in 3D. A projective and adaptively applied pressure estimation is employed to significantly reduce the computational wall time without impacting the overall physical reliability in the modeling of important features of melt segregation, such as melt channel bifurcation in 2D and 3D time dependent simulations.     

Dimension-Reduced Hyperbolic Moment Method for the Boltzmann Equation with BGK-Type Collision

We develop the dimension-reduced hyperbolic moment method for the Boltzmann equation, to improve solution efficiency using a numerical regularized moment method for problems with low-dimensional macroscopic variables and high-dimensional microscopic variables. In the present work, we deduce the globally hyperbolic moment equations for the dimension-reduced Boltzmann equation based on the Hermite expansion and a globally hyperbolic regularization. The numbers of Maxwell boundary condition required for well-posedness are studied. The numerical scheme is then developed and an improved projection algorithm between two different Hermite expansion spaces is developed. By solving several benchmark problems, we validate the method developed and demonstrate the significant efficiency improvement by dimension-reduction.     

A Stochastic Collocation Approach to Bayesian Inference in Inverse Problems

We present an efficient numerical strategy for the Bayesian solution of inverse problems. Stochastic collocation methods, based on generalized polynomial chaos (gPC), are used to construct a polynomial approximation of the forward solution over the support of the prior distribution. This approximation then defines a surrogate posterior probability density that can be evaluated repeatedly at minimal computational cost. The ability to simulate a large number of samples from the posterior distribution results in very accurate estimates of the inverse solution and its associated uncertainty. Combined with high accuracy of the gPC-based forward solver, the new algorithm can provide great efficiency in practical applications. A rigorous error analysis of the algorithm is conducted, where we establish convergence of the approximate posterior to the true posterior and obtain an estimate of the convergence rate. It is proved that fast (exponential) convergence of the gPC forward solution yields similarly fast (exponential) convergence of the posterior. The numerical strategy and the predicted convergence rates are then demonstrated on nonlinear inverse problems of varying smoothness and dimension.     

A Boundary element method based on cauchy integrals for some linear quadratic boundary control problems on a circle

For certain types of elliptic boundary control problems, the boundary element method has considerable advantage over the traditional finite element or finite difference methods because of the reduction of dimensionality in computations. In this paper we examine a variant of such boundary integral methods based on Cauchy integrals. The cost functional here contains only finitely many quadratic terms related to sensory data at those finite interior points. We see that the numerical efficiency of this approach hinges largely on the complexity of the inverse of a certain boundary integral operator. In the case of a circle, such an inverse is readily obtainable and entire computations require only a small effort to yield useful numerical information about the optimal control. Other general situations are also discussed.