Numerical Methods for Solving the Hartree-Fock Equations of Diatomic Molecules I

The theory of domain decomposition is described and used to divide the variable domain of a diatomic molecule into separate regions which are solved independently. This approach makes it possible to use fast Krylov methods in the broad interior of the region while using explicit methods such as Gaussian elimination on the boundaries. As is demonstrated by solving a number of model problems, these methods enable one to obtain solutions of the relevant partial differential equations and eigenvalue equations accurate to six significant figures with a small amount of computational time. Since the numerical approach described in this article decomposes the variable space into separate regions where the equations are solved independently, our approach is very well-suited to parallel computing and offers the long term possibility of studying complex molecules by dividing them into smaller fragments that are calculated separately.     

High Order Finite Difference Discretization for Composite Grid Hierarchy and Its Applications

We introduce efficient approaches to construct high order finite difference discretizations for solving partial differential equations, based on a composite grid hierarchy. We introduce a modification of the traditional point clustering algorithm, obtained by adding restrictive parameters that control the minimal patch length and the size of the buffer zone. As a result, a reduction in the number of interfacial cells is observed. Based on a reasonable geometric grid setting, we discuss a general approach for the construction of stencils in a composite grid environment. The straightforward approach leads to an ill-posed problem. In our approach we regularize this problem, and transform it into solving a symmetric system of linear of equations. Finally, a stencil repository has been designed to further reduce computational overhead. The effectiveness of the discretizations is illustrated by numerical experiments on second order elliptic differential equations.     

A Parallel Domain Decomposition Algorithm for Simulating Blood Flow with Incompressible Navier-Stokes Equations with Resistive Boundary Condition

We introduce and study a parallel domain decomposition algorithm for the simulation of blood flow in compliant arteries using a fully-coupled system of nonlinear partial differential equations consisting of a linear elasticity equation and the incompressible Navier-Stokes equations with a resistive outflow boundary condition. The system is discretized with a finite element method on unstructured moving meshes and solved by a Newton-Krylov algorithm preconditioned with an overlapping restricted additive Schwarz method. The resistive outflow boundary condition plays an interesting role in the accuracy of the blood flow simulation and we provide a numerical comparison of its accuracy with the standard pressure type boundary condition. We also discuss the parallel performance of the implicit domain decomposition method for solving the fully coupled nonlinear system on a supercomputer with a few hundred processors.     

A Tailored Finite Point Method for Solving Steady MHD Duct Flow Problems with Boundary Layers

In this paper we propose a development of the finite difference method, called the tailored finite point method, for solving steady magnetohydrodynamic (MHD) duct flow problems with a high Hartmann number. When the Hartmann number is large, the MHD duct flow is convection-dominated and thus its solution may exhibit localized phenomena such as the boundary layer. Most conventional numerical methods can not efficiently solve the layer problem because they are lacking in either stability or accuracy. However, the proposed tailored finite point method is capable of resolving high gradients near the layer regions without refining the mesh. Firstly, we devise the tailored finite point method for the scalar inhomogeneous convection-diffusion problem, and then extend it to the MHD duct flow which consists of a coupled system of convection-diffusion equations. For each interior grid point of a given rectangular mesh, we construct a finite-point difference operator at that point with some nearby grid points, where the coefficients of the difference operator are tailored to some particular properties of the problem. Numerical examples are provided to show the high performance of the proposed method.     

Extended Physics-Informed Neural Networks (XPINNs): A Generalized Space-Time Domain Decomposition Based Deep Learning Framework for Nonlinear Partial Differential Equations

We propose a generalized space-time domain decomposition approach for the physics-informed neural networks (PINNs) to solve nonlinear partial differential equations (PDEs) on arbitrary complex-geometry domains. The proposed framework, named eXtended PINNs ($XPINNs$), further pushes the boundaries of both PINNs as well as conservative PINNs (cPINNs), which is a recently proposed domain decomposition approach in the PINN framework tailored to conservation laws. Compared to PINN, the XPINN method has large representation and parallelization capacity due to the inherent property of deployment of multiple neural networks in the smaller subdomains. Unlike cPINN, XPINN can be extended to any type of PDEs. Moreover, the domain can be decomposed in any arbitrary way (in space and time), which is not possible in cPINN. Thus, XPINN offers both space and time parallelization, thereby reducing the training cost more effectively. In each subdomain, a separate neural network is employed with optimally selected hyperparameters, e.g., depth/width of the network, number and location of residual points, activation function, optimization method, etc. A deep network can be employed in a subdomain with complex solution, whereas a shallow neural network can be used in a subdomain with relatively simple and smooth solutions. We demonstrate the versatility of XPINN by solving both forward and inverse PDE problems, ranging from one-dimensional to three-dimensional problems, from time-dependent to time-independent problems, and from continuous to discontinuous problems, which clearly shows that the XPINN method is promising in many practical problems. The proposed XPINN method is the generalization of PINN and cPINN methods, both in terms of applicability as well as domain decomposition approach, which efficiently lends itself to parallelized computation. The XPINN code is available on $https://github.com/AmeyaJagtap/XPINNs$.     

An Indirect-Forcing Immersed Boundary Method for Incompressible Viscous Flows with Interfaces on Irregular Domains

An indirect-forcing immersed boundary method for solving the incompressible Navier-Stokes equations involving the interfaces and irregular domains is developed. The rigid boundaries and interfaces are represented by a number of Lagrangian control points. Stationary rigid boundaries are embedded in the Cartesian grid and singular forces at the rigid boundaries are applied to impose the prescribed velocity conditions. The singular forces at the interfaces and the rigid boundaries are then distributed to the nearby Cartesian grid points using the immersed boundary method. In the present work, the singular forces at the rigid boundaries are computed implicitly by solving a small system of equations at each time step to ensure that the prescribed velocity condition at the rigid boundary is satisfied exactly. For deformable interfaces, the forces that the interface exerts on the fluid are computed from the configuration of the elastic interface and are applied to the fluid. The Navier-Stokes equations are discretized using finite difference method on a staggered uniform Cartesian grid by a second order accurate projection method. The ability of the method to simulate viscous flows with interfaces on irregular domains is demonstrated by applying to the rotational flow problem, the relaxation of an elastic membrane and flow in a constriction with an immersed elastic membrane.     

An Improved Integration Scheme for Mode-Coupling-Theory Equations

Within the mode-coupling theory (MCT) of the glass transition, we reconsider the numerical schemes to evaluate the MCT functional. Here we propose nonuniform discretizations of the wave number, in contrast to the standard equidistant grid, in order to decrease the number of grid points without losing accuracy. We discuss in detail how the integration scheme on the new grids has to be modified from standard Riemann integration. We benchmark our approach by solving the MCT equations numerically for mono-disperse hard disks and hard spheres and by computing the critical packing fraction and the nonergodicity parameters. Our results show that significant improvements in performance can be obtained employing a nonuniform grid.     

Numerical Solution of 3D Poisson-Nernst-Planck Equations Coupled with Classical Density Functional Theory for Modeling Ion and Electron Transport in a Confined Environment

We have developed efficient numerical algorithms for solving 3D steady-state Poisson-Nernst-Planck (PNP) equations with excess chemical potentials described by the classical density functional theory (cDFT). The coupled PNP equations are discretized by a finite difference scheme and solved iteratively using the Gummel method with relaxation. The Nernst-Planck equations are transformed into Laplace equations through the Slotboom transformation. Then, the algebraic multigrid method is applied to efficiently solve the Poisson equation and the transformed Nernst-Planck equations. A novel strategy for calculating excess chemical potentials through fast Fourier transforms is proposed, which reduces computational complexity from $\mathcal{O}$($N^2$) to $\mathcal{O}$($NlogN$), where $N$ is the number of grid points. Integrals involving the Dirac delta function are evaluated directly by coordinate transformation, which yields more accurate results compared to applying numerical quadrature to an approximated delta function. Numerical results for ion and electron transport in solid electrolyte for lithium-ion (Li-ion) batteries are shown to be in good agreement with the experimental data and the results from previous studies.     

Explicit Computation of Robin Parameters in Optimized Schwarz Waveform Relaxation Methods for Schrödinger Equations Based on Pseudodifferential Operators

The Optimized Schwarz Waveform Relaxation algorithm, a domain decomposition method based on Robin transmission condition, is becoming a popular computational method for solving evolution partial differential equations in parallel. Along with well-posedness, it offers a good balance between convergence rate, efficient computational complexity and simplicity of the implementation. The fundamental question is the selection of the Robin parameter to optimize the convergence of the algorithm. In this paper, we propose an approach to explicitly estimate the Robin parameter which is based on the approximation of the transmission operators at the subdomain interfaces, for the linear/nonlinear Schrödinger equation. Some illustrating numerical experiments are proposed for the one- and two-dimensional problems.     

A Numerical Comparison Between Quasi-Monte Carlo and Sparse Grid Stochastic Collocation Methods

Quasi-Monte Carlo methods and stochastic collocation methods based on sparse grids have become popular with solving stochastic partial differential equations. These methods use deterministic points for multi-dimensional integration or interpolation without suffering from the curse of dimensionality. It is not evident which method is best, specially on random models of physical phenomena. We numerically study the error of quasi-Monte Carlo and sparse grid methods in the context of ground-water flow in heterogeneous media. In particular, we consider the dependence of the variance error on the stochastic dimension and the number of samples/collocation points for steady flow problems in which the hydraulic conductivity is a lognormal process. The suitability of each technique is identified in terms of computational cost and error tolerance.     

Eulerian Algorithms for Computing the Forward Finite Time Lyapunov Exponent Without Finite Difference upon the Flow Map

In this paper, effective Eulerian algorithms are introduced for the computation of the forward finite time Lyapunov exponent (FTLE) of smooth flow fields. The advantages of the proposed algorithms mainly manifest in two aspects. First, previous Eulerian approaches for computing the FTLE field are improved so that the Jacobian of the flow map can be obtained by directly solving a corresponding system of partial differential equations, rather than by implementing certain finite difference upon the flow map, which can significantly improve the accuracy of the numerical solution especially near the FTLE ridges. Second, in the proposed algorithms, all computations are done on the fly, that is, all required partial differential equations are solved forward in time, which is practically more natural. The new algorithms still maintain the optimal computational complexity as well as the second order accuracy. Numerical examples demonstrate the effectiveness of the proposed algorithms.     

Weak Galerkin Method for Second-Order Elliptic Equations with Newton Boundary Condition

The weak Galerkin (WG) method is a nonconforming numerical method for solving partial differential equations. In this paper, we introduce the WG method for elliptic equations with Newton boundary condition in bounded domains. The Newton boundary condition is a nonlinear boundary condition arising from science and engineering applications. We prove the well-posedness of the WG scheme by the monotone operator theory and the embedding inequality of weak finite element functions. The error estimates are derived. Numerical experiments are presented to verify the theoretical analysis.     

Super-Grid Modeling of the Elastic Wave Equation in Semi-Bounded Domains

We develop a super-grid modeling technique for solving the elastic wave equation in semi-bounded two- and three-dimensional spatial domains. In this method, waves are slowed down and dissipated in sponge layers near the far-field boundaries. Mathematically, this is equivalent to a coordinate mapping that transforms a very large physical domain to a significantly smaller computational domain, where the elastic wave equation is solved numerically on a regular grid. To damp out waves that become poorly resolved because of the coordinate mapping, a high order artificial dissipation operator is added in layers near the boundaries of the computational domain. We prove by energy estimates that the super-grid modeling leads to a stable numerical method with decreasing energy, which is valid for heterogeneous material properties and a free surface boundary condition on one side of the domain. Our spatial discretization is based on a fourth order accurate finite difference method, which satisfies the principle of summation by parts. We show that the discrete energy estimate holds also when a centered finite difference stencil is combined with homogeneous Dirichlet conditions at several ghost points outside of the far-field boundaries. Therefore, the coefficients in the finite difference stencils need only be boundary modified near the free surface. This allows for improved computational efficiency and significant simplifications of the implementation of the proposed method in multi-dimensional domains. Numerical experiments in three space dimensions show that the modeling error from truncating the domain can be made very small by choosing a sufficiently wide super-grid damping layer. The numerical accuracy is first evaluated against analytical solutions of Lamb's problem, where fourth order accuracy is observed with a sixth order artificial dissipation. We then use successive grid refinements to study the numerical accuracy in the more complicated motion due to a point moment tensor source in a regularized layered material.     

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.     

Space-Time Discontinuous Galerkin Method for Maxwell's Equations

A fully discrete discontinuous Galerkin method is introduced for solving time-dependent Maxwell's equations. Distinguished from the Runge-Kutta discontinuous Galerkin method (RKDG) and the finite element time domain method (FETD), in our scheme, discontinuous Galerkin methods are used to discretize not only the spatial domain but also the temporal domain. The proposed numerical scheme is proved to be unconditionally stable, and a convergent rate $\mathcal{O}((∆t)^{r+1}+h^{k+1/2})$ is established under the $L^2$ -norm when polynomials of degree at most $r$ and $k$ are used for temporal and spatial approximation, respectively. Numerical results in both 2-D and 3-D are provided to validate the theoretical prediction. An ultra-convergence of order $(∆t)^{2r+1}$ in time step is observed numerically for the numerical fluxes w.r.t. temporal variable at the grid points.     

Solutions of the 1D Coupled Nonlinear Schrödinger Equations by the CIP-BS Method

In this paper, we present solutions for the one-dimensional coupled nonlinear Schrödinger (CNLS) equations by the Constrained Interpolation Profile-Basis Set (CIP-BS) method. This method uses a simple polynomial basis set, by which physical quantities are approximated with their values and derivatives associated with grid points. Nonlinear operations on functions are carried out in the framework of differential algebra. Then, by introducing scalar products and requiring the residue to be orthogonal to the basis, the linear and nonlinear partial differential equations are reduced to ordinary differential equations for values and spatial derivatives. The method gives stable, less diffusive, and accurate results for the CNLS equations.     

Boundary Control Problems in Convective Heat Transfer with Lifting Function Approach and Multigrid Vanka-Type Solvers

This paper deals with boundary optimal control problems for the heat and Navier-Stokes equations and addresses the issue of defining controls in function spaces which are naturally associated with the volume variables by trace restriction. For this reason we reformulate the boundary optimal control problem into a distributed problem through a lifting function approach. The stronger regularity requirements which are imposed by standard boundary control approaches can then be avoided. Furthermore, we propose a new numerical strategy that allows solving the coupled optimality system in a robust way for a large number of unknowns. The optimality system resulting from a finite element discretization is solved by a local multigrid algorithm with domain decomposition Vanka-type smoothers. The purpose of these smoothers is to solve the optimality system implicitly over subdomains with a small number of degrees of freedom, in order to achieve robustness with respect to the regularization parameters in the cost functional. We present the results of some test cases where temperature is the observed quantity and the control quantity corresponds to the boundary values of the fluid temperature in a portion of the boundary. The control region for the observed quantity is a part of the domain where it is interesting to match a desired temperature value.     

An Augmented Approach for the Pressure Boundary Condition in a Stokes Flow

An augmented method is proposed for solving stationary incompressible Stokes equations with a Dirichlet boundary condition along parts of the boundary. In this approach, the normal derivative of the pressure along the parts of the boundary is introduced as an additional variable and it is solved by the GMRES iterative method. The dimension of the augmented variable in discretization is the number of grid points along the boundary which is O(N). Each GMRES iteration (or one matrix-vector multiplication) requires three fast Poisson solvers for the pressure and the velocity. In our numerical experiments, only a few iterations are needed. We have also combined the augmented approach for Stokes equations involving interfaces, discontinuities, and singularities.     

An Efficient,Energy Stable Scheme for the Cahn-Hilliard-Brinkman System

We present an unconditionally energy stable and uniquely solvable finite difference scheme for the Cahn-Hilliard-Brinkman (CHB) system, which is comprised of a Cahn-Hilliard-type diffusion equation and a generalized Brinkman equation modeling fluid flow. The CHB system is a generalization of the Cahn-Hilliard-Stokes model and describes two phase very viscous flows in porous media. The scheme is based on a convex splitting of the discrete CH energy and is semi-implicit. The equations at the implicit time level are nonlinear, but we prove that they represent the gradient of a strictly convex functional and are therefore uniquely solvable, regardless of time step size. Owing to energy stability, we show that the scheme is stable in the time and space discrete$ℓ^∞$(0,$T$;$H^1_h$) and $ℓ^2$(0,$T$;$H^2_h$) norms. We also present an efficient, practical nonlinear multigrid method – comprised of a standard FAS method for the Cahn-Hilliard part, and a method based on the Vanka smoothing strategy for the Brinkman part – for solving these equations. In particular, we provide evidence that the solver has nearly optimal complexity in typical situations. The solver is applied to simulate spinodal decomposition of a viscous fluid in a porous medium, as well as to the more general problems of buoyancy- and boundary-driven flows.     

舟月骨间韧带在腕尺桡偏时舟月骨三维运动中的作用

目的 观察在腕关节尺桡偏运动时,舟月骨间韧带(scapholunate interosseous ligament,SLIL)及其各亚区分级切断后对舟骨与月骨三维运动的影响,从亚区水平对SLIL在舟月骨三维运动中的作用进行研究.方法 取12例新鲜冰冻成人上肢标本(左、右侧各6例),共分为正常对照组、断SLIL近侧组、断SLIL近+背侧组、断SLIL近+掌侧组、SLIL全断组(每组6例).采用三维激光扫描及图像重建技术,测定在腕关节尺桡偏时舟骨与月骨的三维运动数据,并进行统计学处理.结果 在正常腕关节做尺桡偏运动时,舟、月骨同时也做尺桡偏运动,但还伴有明显的屈伸方向上的运动;与此同时,舟、月骨还做一定的旋前、旋后活动.在切断近+掌和近+背侧亚区组,舟月骨的活动受到影响.全切断SLIL会导致舟骨掌屈度增加,月骨掌屈活动减弱.结论 采用三维激光扫描及图像重建技术可准确测定腕舟月骨的三维活动度,从SLIL亚区水平上进行的研究表明,部分或全部切断SLIL后会对舟月骨的掌屈程度产生明显影响,但SLIL近侧亚区在控制舟、月骨的活动上不起主要作用.