Numerical oscillation of the total energy can be observed when the Kohn- Sham equation is solved by real-space methods to simulate the translational move of an electronic system. Effectively remove or reduce the unphysical oscillation is crucial not only for the optimization of the geometry of the electronic structure, but also for the study of molecular dynamics. In this paper, we study such unphysical oscillation based on the numerical framework in [G. Bao, G. H. Hu, and D. Liu, An h-adaptive finite element solver for the calculations of the electronic structures, Journal of Computational Physics, Volume 231, Issue 14, Pages 4967–4979, 2012], and deliver some numerical methods to constrain such unphysical effect for both pseudopotential and all-electron calculations, including a stabilized cubature strategy for Hamiltonian operator, and an a posteriori error estimator of the finite element methods for Kohn-Sham equation. The numerical results demonstrate the effectiveness of our method on restraining unphysical oscillation of the total energies.

In this paper, we study the Camassa-Holm equation and the Degasperis-Procesi equation. The two equations are in the family of integrable peakon equations, and both have very rich geometric properties. Based on these geometric structures, we construct the geometric numerical integrators for simulating their soliton solutions. The Camassa-Holm equation and the Degasperis-Procesi equation have many common properties, however they also have the significant difference, for example there exist the shock wave solutions for the Degasperis-Procesi equation. By using the symplectic Fourier pseudo-spectral integrator, we simulate the peakon solutions of the two equations. To illustrate the smooth solitons and shock wave solutions of the DP equation, we use the splitting technique and combine the composition methods. In the numerical experiments, comparisons of these two kinds of methods are presented in terms of accuracy, computational cost and invariants preservation.

Transient diffusion equations arise in many branches of engineering and applied sciences (e.g., heat transfer and mass transfer), and are parabolic partial differential equations. It is well-known that these equations satisfy important mathematical properties like maximum principles and the non-negative constraint, which have implications in mathematical modeling. However, existing numerical formulations for these types of equations do not, in general, satisfy maximum principles and the non-negative constraint. In this paper, we present a methodology for enforcing maximum principles and the non-negative constraint for transient anisotropic diffusion equation. The proposed methodology is based on the method of horizontal lines in which the time is discretized first. This results in solving steady anisotropic diffusion equation with decay equation at every discrete time-level. We also present other plausible temporal discretizations, and illustrate their shortcomings in meeting maximum principles and the non-negative constraint. The proposed methodology can handle general computational grids with no additional restrictions on the time-step. We illustrate the performance and accuracy of the proposed methodology using representative numerical examples. We also perform a numerical convergence analysis of the proposed methodology. For comparison, we also present the results from the standard single-field semi-discrete formulation and the results from a popular software package, which all will violate maximum principles and the non-negative constraint.

An efficient and accurate numerical scheme is proposed for solving the transverse electric (TE) mode electromagnetic (EM) propagation problem in two-dimensional earth. The scheme is based on the alternating direction finite-difference time-domain (ADI-FDTD) method. Unlike the conventional upward continuation approach for the earth-air interface, an integral formulation for the interface boundary is developed and it is effectively incorporated to the ADI solver. Stability and convergence analysis together with an error estimate are presented. Numerical simulations are carried out to validate the proposed method, and the advantage of the present method over the popular Du-Fort-Frankel scheme is clearly demonstrated. Examples of the electromagnetic field propagation in the ground with anomaly further verify the effectiveness of the proposed scheme.

A front tracking method combined with the real ghost fluid method (RGFM) is proposed for simulations of fluid interfaces in two-dimensional compressible flows. In this paper the Riemann problem is constructed along the normal direction of interface and the corresponding Riemann solutions are used to track fluid interfaces. The interface boundary conditions are defined by the RGFM, and the fluid interfaces are explicitly tracked by several connected marker points. The Riemann solutions are also used directly to update the flow states on both sides of the interface in the RGFM. In order to validate the accuracy and capacity of the new method, extensive numerical tests including the bubble advection, the Sod tube, the shock-bubble interaction, the Richtmyer-Meshkov instability and the gas-water interface, are simulated by using the Euler equations. The computational results are also compared with earlier computational studies and it shows good agreements including the compressible gas-water system with large density differences.

In this paper we consider PDE-constrained optimization problems which incorporate an H1 regularization control term. We focus on a time-dependent PDE, and consider both distributed and boundary control. The problems we consider include bound constraints on the state, and we use a Moreau-Yosida penalty function to handle this. We propose Krylov solvers and Schur complement preconditioning strategies for the different problems and illustrate their performance with numerical examples.

Because of stability constraints, most numerical schemes applied to hyperbolic systems of equations turn out to be costly when the flux term is multiplied by some very large scalar. This problem emerges with the M1 system of equations in the field of radiotherapy when considering heterogeneous media with very disparate densities. Additionally, the flux term of the M1 system is non-linear, and in order for the model to be well-posed the numerical solution needs to fulfill conditions called realizability. In this paper, we propose a numerical method that overcomes the stability constraint and preserves the realizability property. For this purpose, we relax the M1 system to obtain a linear flux term. Then we extend the stencil of the difference quotient to obtain stability. The scheme is applied to a radiotherapy dose calculation example.

We devise an efficient algorithm for the symbolic calculation of irreducible angular momentum and spin (LS) eigenspaces within the n-fold antisymmetrized tensor product ^n Vu, where n is the number of electrons and u = s,p,d,… denotes the atomic subshell. This is an essential step for dimension reduction in configuration-interaction (CI) methods applied to atomic many-electron quantum systems. The algorithm relies on the observation that each Lz eigenstate with maximal eigenvalue is also an L2 eigenstate (equivalently for Sz and S2), as well as the traversal of LS eigenstates using the lowering operators L_ and S_. Iterative application to the remaining states in ^nVu leads to an implicit simultaneous diagonalization. A detailed complexity analysis for fixed n and increasing subshell number u yields run time O(u3n–2). A symbolic computer algebra implementation is available online.

The massively parallel, nonlinear, three-dimensional (3D), toroidal, electrostatic, gyrokinetic, particle-in-cell (PIC), Cartesian geometry UCAN code, with particle ions and adiabatic electrons, has been successfully exercised to identify non-diffusive transport characteristics in present day tokamak discharges. The limitation in applying UCAN to larger scale discharges is the 1D domain decomposition in the toroidal (or z-) direction for massively parallel implementation using MPI which has restricted the calculations to a few hundred ion Larmor radii or gyroradii per plasma minor radius. To exceed these sizes, we have implemented 2D domain decomposition in UCAN with the addition of the y-direction to the processor mix. This has been facilitated by use of relevant components in the P2LIB library of field and particle management routines developed for UCLA's UPIC Framework of conventional PIC codes. The gyro-averaging specific to gyrokinetic codes is simplified by the use of replicated arrays for efficient charge accumulation and force deposition. The 2D domain-decomposed UCAN2 code reproduces the original 1D domain nonlinear results within round-off. Benchmarks of UCAN2 on the Cray XC30 Edison at NERSC demonstrate ideal scaling when problem size is increased along with processor number up to the largest power of 2 available, namely 131,072 processors. These particle weak scaling benchmarks also indicate that the 1 nanosecond per particle per time step and 1 TFlops barriers are easily broken by UCAN2 with 1 billion particles or more and 2000 or more processors.

In this paper, we introduce and study a new method for solving inverse source problems, through a working model that arises in bioluminescence tomography (BLT). In the BLT problem, one constructs quantitatively the bioluminescence source distribution inside a small animal from optical signals detected on the animal's body surface. The BLT problem possesses strong ill-posedness and often the Tikhonov regularization is used to obtain stable approximate solutions. In conventional Tikhonov regularization, it is crucial to choose a proper regularization parameter for trade off between the accuracy and stability of approximate solutions. The new method is based on a combination of the boundary condition and the boundary measurement in a parameter-dependent single complex Robin boundary condition, followed by the Tikhonov regularization. By properly adjusting the parameter in the Robin boundary condition, we achieve two important properties for our new method: first, the regularized solutions are uniformly stable with respect to the regularization parameter so that the regularization parameter can be chosen based solely on the consideration of the solution accuracy; second, the convergence order of the regularized solutions reaches one with respect to the noise level. Then, the finite element method is used to compute numerical solutions and a new finite element error estimate is derived for discrete solutions. These results improve related results found in the existing literature. Several numerical examples are provided to illustrate the theoretical results.

Using the incompressible isotropic turbulent fields obtained from direct numerical simulation and large-eddy simulation, we studied the statistics of oscillation structures based on local zero-crossings and their relation with inertial-range intermittency for transverse velocity and passive scalar. Our results show that for both the velocity and passive scalar, the local oscillation structures are statistically scale-invariant at high Reynolds number, and the inertial-range intermittency of the overall flow region is determined by the most intermittent structures characterized by one local zero-crossing. Local flow patterns conditioned on the oscillation structures are characterized by the joint probability density function of the invariants of the filtered velocity gradient tensor at inertial range. We demonstrate that the most intermittent regions for longitudinal velocity tend to lay at the saddle area, while those for the transverse velocity tend to locate at the vortex-dominated area. The connection between the ramp-cliff structures in passive scalar field and the corresponding saddle regions in the velocity field is also verified by the approach of oscillation structure classification.