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.     

A Third Order BDF Energy Stable Linear Scheme for the No-Slope-Selection Thin Film Model

In this paper we propose and analyze a (temporally) third order accurate backward differentiation formula (BDF) numerical scheme for the no-slope-selection (NSS) equation of the epitaxial thin film growth model, with Fourier pseudo-spectral discretization in space. The surface diffusion term is treated implicitly, while the nonlinear chemical potential is approximated by a third order explicit extrapolation formula for the sake of solvability. In addition, a third order accurate Douglas-Dupont regularization term, in the form of $−A∆t^2∆^2_N (u^{n+1}−u^n)$, is added in the numerical scheme. A careful energy stability estimate, combined with Fourier eigenvalue analysis, results in the energy stability in a modified version, and a theoretical justification of the coefficient $A$ becomes available. As a result of this energy stability analysis, a uniform in time bound of the numerical energy is obtained. And also, the optimal rate convergence analysis and error estimate are derived in details, in the $ℓ^∞(0,T;ℓ^2)∩ℓ^2(0,T;H^2_h)$ norm, with the help of a linearized estimate for the nonlinear error terms. Some numerical simulation results are presented to demonstrate the efficiency of the numerical scheme and the third order convergence. The long time simulation results for $ε = 0.02$ (up to $T = 3×10^5$) have indicated a logarithm law for the energy decay, as well as the power laws for growth of the surface roughness and the mound width. In particular, the power index for the surface roughness and the mound width growth, created by the third order numerical scheme, is more accurate than those produced by certain second order energy stable schemes in the existing literature.     

A Modified Crank-Nicolson Numerical Scheme for the Flory-Huggins Cahn-Hilliard Model

In this paper we propose and analyze a second order accurate numerical scheme for the Cahn-Hilliard equation with logarithmic Flory Huggins energy potential. A modified Crank-Nicolson approximation is applied to the logarithmic nonlinear term, while the expansive term is updated by an explicit second order Adams-Bashforth extrapolation, and an alternate temporal stencil is used for the surface diffusion term. A nonlinear artificial regularization term is added in the numerical scheme, which ensures the positivity-preserving property, i.e., the numerical value of the phase variable is always between -1 and 1 at a point-wise level. Furthermore, an unconditional energy stability of the numerical scheme is derived, leveraging the special form of the logarithmic approximation term. In addition, an optimal rate convergence estimate is provided for the proposed numerical scheme, with the help of linearized stability analysis. A few numerical results, including both the constant-mobility and solution-dependent mobility flows, are presented to validate the robustness of the proposed numerical scheme.     

An Adaptive Finite Element Method with Hybrid Basis for Singularly Perturbed Nonlinear Eigenvalue Problems

In this paper, we propose a uniformly convergent adaptive finite element method with hybrid basis (AFEM-HB) for the discretization of singularly perturbed nonlinear eigenvalue problems under constraints with applications in Bose-Einstein condensation (BEC) and quantum chemistry. We begin with the time-independent Gross-Pitaevskii equation and show how to reformulate it into a singularly perturbed nonlinear eigenvalue problem under a constraint. Matched asymptotic approximations for the problem are reviewed to confirm the asymptotic behaviors of the solutions in the boundary/interior layer regions. By using the normalized gradient flow, we propose an adaptive finite element with hybrid basis to solve the singularly perturbed nonlinear eigenvalue problem. Our basis functions and the mesh are chosen adaptively to the small parameter ε. Extensive numerical results are reported to show the uniform convergence property of our method. We also apply the AFEM-HB to compute the ground and excited states of BEC with box/harmonic/optical lattice potential in the semiclassical regime (0<ε≪1). In addition, we give a detailed error analysis of our AFEM-HB to a simpler singularly perturbed two point boundary value problem, show that our method has a minimum uniform convergence order $\mathcal{O}$(1/$(NlnN)^\frac{2}{3}$).     

A Frequency Determination Method for Digitized NMR Signals

We present a high precision frequency determination method for digitized NMR FID signals. The method employs high precision numerical integration rather than simple summation as in many other techniques. With no independent knowledge of the other parameters of a NMR FID signal (phase $ϕ$, amplitude $A$, and transverse relaxation time $T_2$) this method can determine the signal frequency $f_0$ with a precision of 1/(8$π^2$$f^2_0$$T^2_2$) if the observation time $T$ ≫$T_2$. The method is especially convenient when the detailed shape of the observed FT NMR spectrum is not well defined. When $T_2$ is +∞ and the signal becomes pure sinusoidal, the precision of the method is 3/(2$π^2$$f^2_0$$T_2$) which is one order more precise than the ±1 count error induced precision of a typical frequency counter. Analysis of this method shows that the integration reduces the noise by bandwidth narrowing as in a lock-in amplifier, and no extra signal filters are needed. For a pure sinusoidal signal we find from numerical simulations that the noise-induced error in this method reaches the Cramer-Rao Lower Band (CRLB) on frequency determination. For the damped sinusoidal case of most interest, the noise-induced error is found to be within a factor of 2 of CRLB when the measurement time $T$ is 2 or 3 times larger than $T_2$. We discuss possible improvements for the precision of this method.     

Decoupled,Energy Stable Numerical Scheme for the Cahn-Hilliard-Hele-Shaw System with Logarithmic Flory-Huggins Potential

In this paper, a decoupling numerical method for solving Cahn-Hilliard-Hele-Shaw system with logarithmic potential is proposed. Combing with a convex-splitting of the energy functional, the discretization of the Cahn-Hilliard equation in time is presented. The nonlinear term in Cahn-Hilliard equation is decoupled from the pressure gradient by using a fractional step method. Therefore, to update the pressure, we just need to solve a Possion equation at each time step by using an incremental pressure-correction technique for the pressure gradient in Darcy equation. For logarithmic potential, we use the regularization procedure, which make the domain for the regularized functional $F$($ф$) is extended from (−1,1) to (−∞,∞). Further, the stability and the error estimate of the proposed method are proved. Finally, a series of numerical experiments are implemented to illustrate the theoretical analysis.     

A Decreasing Upper Bound of the Energy for Time-Fractional Phase-Field Equations

In this article, we study the energy dissipation property of time-fractional Allen–Cahn equation. On the continuous level, we propose an upper bound of energy that decreases with respect to time and coincides with the original energy at $t=0$ and as $t$ tends to $ ∞.$ This upper bound can also be viewed as a nonlocal-in-time modified energy which is the summation of the original energy and an accumulation term due to the memory effect of time-fractional derivative. In particular, the decrease of the modified energy indicates that the original energy indeed decays w.r.t. time in a small neighborhood at $t=0.$ We illustrate the theory mainly with the time-fractional Allen–Cahn equation but it could also be applied to other time-fractional phase-field models such as the Cahn–Hilliard equation. On the discrete level, the decreasing upper bound of energy is useful for proving energy dissipation of numerical schemes. First-order L1 and second-order L2 schemes for the time-fractional Allen–Cahn equation have similar decreasing modified energies, so that stability can be established. Some numerical results are provided to illustrate the behavior of this modified energy and to verify our theoretical results.     

A Mixed Finite Element Scheme for Biharmonic Equation with Variable Coefficient and von Kármán Equations

In this paper, a new mixed finite element scheme using element-wise stabilization is introduced for the biharmonic equation with variable coefficient on Lipschitz polyhedral domains. The proposed scheme doesn’t involve any integration along mesh interfaces. The gradient of the solution is approximated by $H$(div)-conforming $BDM_{k+1}$ element or vector valued Lagrange element with order $k+1,$ while the solution is approximated by Lagrange element with order $k+2$ for any $k≥0.$ This scheme can be easily implemented and produces symmetric and positive definite linear system. We provide a new discrete $H^2$-norm stability, which is useful not only in analysis of this scheme but also in ${\rm C}^0$ interior penalty methods and DG methods. Optimal convergences in both discrete $H^2$-norm and $L^2$-norm are derived. This scheme with its analysis is further generalized to the von Kármán equations. Finally, numerical results verifying the theoretical estimates of the proposed algorithms are also presented.     

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.     

An All Speed Second Order IMEX Relaxation Scheme for the Euler Equations

We present an implicit-explicit finite volume scheme for the Euler equations. We start from the non-dimensionalised Euler equations where we split the pressure in a slow and a fast acoustic part. We use a Suliciu type relaxation model which we split in an explicit part, solved using a Godunov-type scheme based on an approximate Riemann solver, and an implicit part where we solve an elliptic equation for the fast pressure. The relaxation source terms are treated projecting the solution on the equilibrium manifold. The proposed scheme is positivity preserving with respect to the density and internal energy and asymptotic preserving towards the incompressible Euler equations. For this first order scheme we give a second order extension which maintains the positivity property. We perform numerical experiments in 1D and 2D to show the applicability of the proposed splitting and give convergence results for the second order extension.     

Exponential Runge-Kutta Methods for the Multispecies Boltzmann Equation

This paper generalizes the exponential Runge-Kutta asymptotic preserving (AP) method developed in [G. Dimarco and L. Pareschi, SIAM Numer. Anal., 49 (2011), pp. 2057–2077] to compute the multi-species Boltzmann equation. Compared to the single species Boltzmann equation that the method was originally applied to, this set of equation presents a new difficulty that comes from the lack of local conservation laws due to the interaction between different species. Hence extra stiff nonlinear source terms need to be treated properly to maintain the accuracy and the AP property. The method we propose does not contain any nonlinear nonlocal implicit solver, and can capture the hydrodynamic limit with time step and mesh size independent of the Knudsen number. We prove the positivity and strong AP properties of the scheme, which are verified by two numerical examples.     

A Fast Time Splitting Finite Difference Approach to Gross–Pitaevskii Equations

We propose an idea to solve the Gross–Pitaevskii equation for dark structures inside an infinite constant background density $ρ_∞$=${|ψ_∞|}^2$, without the introduction of artificial boundary conditions. We map the unbounded physical domain $\mathbb{R}^3$ into the bounded domain ${(−1,1)}^3$ and discretize the rescaled equation by equispaced 4th-order finite differences. This results in a free boundary approach, which can be solved in time by the Strang splitting method. The linear part is solved by a new, fast approximation of the action of the matrix exponential at machine precision accuracy, while the nonlinear part can be solved exactly. Numerical results confirm existing ones based on the Fourier pseudospectral method and point out some weaknesses of the latter such as the need of a quite large computational domain, and thus a consequent critical computational effort, in order to provide reliable time evolution of the vortical structures, of their reconnections, and of integral quantities like mass, energy, and momentum. The free boundary approach reproduces them correctly, also in finite subdomains, at low computational cost. We show the versatility of this method by carrying out one- and three-dimensional simulations and by using it also in the case of Bose–Einstein condensates, for which $ψ$→0 as the spatial variables tend to infinity.     

Optimal Error Estimates of Compact Finite Difference Discretizations for the Schrödinger-Poisson System

We study compact finite difference methods for the Schrödinger-Poisson equation in a bounded domain and establish their optimal error estimates under proper regularity assumptions on wave function $ψ$ and external potential $V(x)$. The Crank-Nicolson compact finite difference method and the semi-implicit compact finite difference method are both of order $\mathcal{O}$($h^4$+$τ^2$) in discrete $l^2$, $H^1$ and $l^∞$ norms with mesh size $h$ and time step $τ$. For the errors of compact finite difference approximation to the second derivative and Poisson potential are nonlocal, thus besides the standard energy method and mathematical induction method, the key technique in analysis is to estimate the nonlocal approximation errors in discrete $l^∞$ and $H^1$ norm by discrete maximum principle of elliptic equation and properties of some related matrix. Also some useful inequalities are established in this paper. Finally, extensive numerical results are reported to support our error estimates of the numerical methods.     

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.     

Natural Convection Heat Transfer in a Porous Cavity with Sinusoidal Temperature Distribution Using Cu/Water Nanofluid: Double MRT Lattice Boltzmann Method

In this study, natural convection flow in a porous cavity with sinusoidal temperature distribution has been analyzed by a new double multi relaxation time (MRT) Lattice Boltzmann method (LBM). We consider a copper/water nanofluid filling a porous cavity. For simulating the temperature and flow fields, D2Q5 and D2Q9 lattices are utilized respectively, and the effects of different Darcy numbers (Da) (0.001-0.1) and various Rayleigh numbers (Ra) ($10^3$-$10^5$) for porosity ($ε$) between 0.4 and 0.9 have been considered. Phase deviation ($θ$) changed from 0 to $π$ and the volume fraction of nanoparticles (Ø) varied from 0 to 6%. The present results show a good agreement with the previous works, thus confirming the reliability the new numerical method proposed in this paper. It is indicated that the heat transfer rate increases at increasing Darcy number, porosity, Rayleigh number, the volume fraction of nanoparticles and phase deviation. However, the most sensitive parameter is the Rayleigh number. The maximum Nusselt deviation is 10%, 32% and 33% for Ra=$10^3$, $10^4$ and $10^5$, respectively, with $ε = 0.4$ to $ε = 0.9$. It can be concluded that the effect of Darcy number on the heat transfer rate increases at increasing Rayleigh number, yielding a maximum enhancement of the average Nusselt number around 12% and 61% for Ra=$10^3$ and Ra=$10^5$, respectively.     

Fully Decoupled,Linear and Unconditionally Energy Stable Schemes for the Binary Fluid-Surfactant Model

Here, we develop a first and a second order time stepping schemes for a binary fluid-surfactant phase field model by using the scalar auxiliary variable approach. The free energy contains a double-well potential, a nonlinear coupling entropy and a Flory-Huggins potential. The resulting coupled system consists of a Cahn-Hilliard type equation and a Wasserstein type equation which leads to a degenerate problem. By introducing only one scalar auxiliary variable, the system is transformed into an equivalent form so that the nonlinear terms can be treated semi-explicitly. Both the schemes are linear and decoupled, thus they can be solved efficiently. We further prove that these semi-discretized schemes in time are unconditionally energy stable. Some numerical experiments are performed to validate the accuracy and energy stability of the proposed schemes.     

On Invariant-Preserving Finite Difference Schemes for the Camassa-Holm Equation and the Two-Component Camassa-Holm System

The purpose of this paper is to develop and test novel invariant-preserving finite difference schemes for both the Camassa-Holm (CH) equation and one of its 2-component generalizations (2CH). The considered PDEs are strongly nonlinear, admitting soliton-like peakon solutions which are characterized by a slope discontinuity at the peak in the wave shape, and therefore suitable for modeling both short wave breaking and long wave propagation phenomena. The proposed numerical schemes are shown to preserve two invariants, momentum and energy, hence numerically producing wave solutions with smaller phase error over a long time period than those generated by other conventional methods. We first apply the scheme to the CH equation and showcase the merits of considering such a scheme under a wide class of initial data. We then generalize this scheme to the 2CH equation and test this scheme under several types of initial data.     

Convergence Rate Analysis for Deep Ritz Method

Using deep neural networks to solve PDEs has attracted a lot of attentions recently. However, why the deep learning method works is falling far behind its empirical success. In this paper, we provide a rigorous numerical analysis on deep Ritz method (DRM) [47] for second order elliptic equations with Neumann boundary conditions. We establish the first nonasymptotic convergence rate in $H^1$ norm for DRM using deep networks with ${\rm ReLU}^2$ activation functions. In addition to providing a theoretical justification of DRM, our study also shed light on how to set the hyperparameter of depth and width to achieve the desired convergence rate in terms of number of training samples. Technically, we derive bound on the approximation error of deep ${\rm ReLU}^2$ network in $C^1$ norm and bound on the Rademacher complexity of the non-Lipschitz composition of gradient norm and ${\rm ReLU}^2$ network, both of which are of independent interest.     

Directional $\mathcal{H}^2$ Compression Algorithm: Optimisations and Application to a Discontinuous Galerkin BEM for the Helmholtz Equation

This study aimed to specialise a directional $\mathcal{H}^2 (\mathcal{D}\mathcal{H}^2)$ compression to matrices arising from the discontinuous Galerkin (DG) discretisation of the hypersingular equation in acoustics. The significant finding is an algorithm that takes a DG stiffness matrix and finds a near-optimal $\mathcal{D}\mathcal{H}^2$ approximation for low and high-frequency problems. We introduced the necessary special optimisations to make this algorithm more efficient in the case of a DG stiffness matrix. Moreover, an automatic parameter tuning strategy makes it easy to use and versatile. Numerical comparisons with a classical Boundary Element Method (BEM) show that a DG scheme combined with a $\mathcal{D}\mathcal{H}^2$ gives better computational efficiency than a classical BEM in the case of high-order finite elements and $hp$ heterogeneous meshes. The results indicate that DG is suitable for an auto-adaptive context in integral equations.     

Reliability Investigation of BiCGStab and IDR Solvers for the Advection-Diffusion-Reaction Equation

The reliability of BiCGStab and IDR solvers for the exponential scheme discretization of the advection-diffusion-reaction equation is investigated. The resulting discretization matrices have real eigenvalues. We consider BiCGStab, IDR$(S),$ BiCGStab$(L)$ and various modifications of BiCGStab, where $S$ denotes the dimension of the shadow space and $L$ the degree of the polynomial used in the polynomial part. Several implementations of BiCGStab exist which are equivalent in exact arithmetic, however, not in finite precision arithmetic. The modifications of BiCGStab we consider are; choosing a random shadow vector, a reliable updating scheme, and storing the best intermediate solution. It is shown that the Local Minimal Residual algorithm, a method similar to the "minimize residual" step of BiCGStab, can be interpreted in terms of a time-dependent advection-diffusion-reaction equation with homogeneous Dirichlet boundary conditions for the residual, which plays a key role in the convergence analysis. Due to the real eigenvalues, the benefit of BiCGStab$(L)$ compared to BiCGStab is shown to be modest in numerical experiments. Non-sparse (e.g. uniform random) shadow residual turns out to be essential for the reliability of BiCGStab. The reliable updating scheme ensures the required tolerance is truly achieved. Keeping the best intermediate solution has no significant effect. Recommendation is to modify BiCGStab with a random shadow residual and the reliable updating scheme, especially in the regime of large Péclet and small Damköhler numbers. An alternative option is IDR($S$), which outperforms BiCGStab for problems with strong advection in terms of the number of matrix-vector products. The MATLAB code used in the numerical experiments is available on GitLab: https://gitlab.com/ChrisSchoutrop/krylov-adr, a C++ implementation of IDR$(S)$ is available in the Eigen linear algebra library: http: //eigen.tuxfamily.org.