The local one-dimensional multisymplectic scheme (LOD-MS) is developed for the three-dimensional (3D) Gross-Pitaevskii (GP) equation in Bose-Einstein condensates. The idea is originated from the advantages of multisymplectic integrators and from the cheap computational cost of the local one-dimensional (LOD) method. The 3D GP equation is split into three linear LOD Schrödinger equations and an exactly solvable nonlinear Hamiltonian ODE. The three linear LOD Schrödinger equations are multisymplectic which can be approximated by multisymplectic integrator (MI). The conservative properties of the proposed scheme are investigated. It is mass-preserving. Surprisingly, the scheme preserves the discrete local energy conservation laws and global energy conservation law if the wave function is variable separable. This is impossible for conventional MIs in nonlinear Hamiltonian context. The numerical results show that the LOD-MS can simulate the original problems very well. They are consistent with the numerical analysis.

In this article we detail the methodology developed to construct an efficient interface description technique — the robust conservative level set (RCLS) — to simulate multiphase flows on mixed-element unstructured meshes while conserving mass to machine accuracy. The approach is tailored specifically for industry as the three-dimensional unstructured approach allows for the treatment of very complex geometries. In addition, special care has been taken to optimise the trade-off between accuracy and computational cost while maintaining the robustness of the numerical method. This was achieved by solving the transport equations for the liquid volume fraction using a WENO scheme for polyhedral meshes and by adding a flux-limiter algorithm. The performance of the resulting method has been compared against established multiphase numerical methods and its ability to capture the physics of multiphase flows is demonstrated on a range of relevant test cases. Finally, the RCLS method has been applied to the simulation of the primary breakup of a flat liquid sheet of kerosene in co-flowing high-pressure gas. This quasi-DNS/LES computation was performed at relevant aero-engine conditions on a three-dimensional mixed-element unstructured mesh. The numerical results have been validated qualitatively against theoretical predictions and experimental data. In particular, the expected breakup regime was observed in the simulation results. Finally, the computation reproduced faithfully the breakup length predicted by a correlation based on experimental data. This constitutes a first step towards a quantitative validation.

We propose a new Absorbing Boundary Condition (ABC) for the acoustic wave equation which is derived from a micro-local diagonalization process formerly defined by M.E. Taylor and which does not depend on the geometry of the surface bearing the ABC. By considering the principal symbol of the wave equation both in the hyperbolic and the elliptic regions, we show that a second-order ABC can be constructed as the combination of an existing first-order ABC and a Fourier-Robin condition. We compare the new ABC with other ABCs and we show that it performs well in simple configurations and that it improves the accuracy of the numerical solution without increasing the computational burden.

The gradient of the fluid stresses exerted on curved boundaries, conventionally computed in terms of directional derivatives of a tensor, is here analyzed by using the notion of intrinsic derivative which represents the geometrically appropriate tool for measuring tensor variations projected on curved surfaces. Relevant differences in the two approaches are found by using the classical Stokes analytical solution for the slow motion of a fluid over a fixed sphere and a numerically generated three dimensional dynamical scenario. Implications for theoretical fluid dynamics and for applied sciences are finally discussed.

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.

In this paper, an immersed interface method is presented to simulate the dynamics of inextensible interfaces in an incompressible flow. The tension is introduced as an augmented variable to satisfy the constraint of interface inextensibility and the resulting augmented system is solved by the GMRES method. In this work, the arclength of the interface is locally and globally conserved as the enclosed region undergoes deformation. The forces at the interface are calculated from the configuration of the interface and the computed augmented variable, and then applied to the fluid through the related jump conditions. The governing equations are discretized on a MAC grid via a second-order finite difference scheme which incorporates jump contributions and solved by the conjugate gradient Uzawa-type method. The proposed method is applied to several examples including the deformation of a liquid capsule with inextensible interfaces in a shear flow. Numerical results reveal that both the area enclosed by interface and arclength of interface are conserved well simultaneously. These provide further evidence on the capability of the present method to simulate incompressible flows involving inextensible interfaces.

We extend the full wavefield modeling with forward scattering theory and Volterra Renormalization to a vertically varying two-parameter (velocity and density) acoustic medium. The forward scattering series, derived by applying Born-Neumann iterative procedure to the Lippmann-Schwinger equation (LSE), is a well known tool for modeling and imaging. However, it has limited convergence properties depending on the strength of contrast between the actual and reference medium or the angle of incidence of a plane wave component. Here, we introduce the Volterra renormalization technique to the LSE. The renormalized LSE and related Neumann series are absolutely convergent for any strength of perturbation and any incidence angle. The renormalized LSE can further be separated into two sub-Volterra type integral equations, which are then solved noniteratively. We apply the approach to velocity-only, density-only, and both velocity and density perturbations. We demonstrate that this Volterra Renormalization modeling is a promising and efficient method. In addition, it can also provide insight for developing a scattering theory-based direct inversion method.

In this work, we introduce an IMEX discontinuous Galerkin solver for the weakly compressible isentropic Euler equations. The splitting needed for the IMEX temporal integration is based on the recently introduced reference solution splitting [32, 52], which makes use of the incompressible solution. We show that the overall method is asymptotic preserving. Numerical results show the performance of the algorithm which is stable under a convective CFL condition and does not show any order degradation.

An alternate BGK type formulation of the Enskog equation has been recently proposed. It was shown that the new model has a valid H-theorem and correct thermal conductivity. We propose Lattice Boltzmann (LB) formulation of this new Enskog-BGK model. The molecular nature of the model is verified in case of shear flow by comparing the predicted normal stress behavior by the current model with the prediction of molecular dynamics simulations. We extend the model for multiphase flow by incorporating attractive part as Vlasov type force. To validate multiphase formulation, the results of 3D simulations of a condensing bubble in a periodic box are presented.

We propose a mathematically rigorous method to measure the spontaneous curvature of a bilayer membrane by molecular dynamics (MD) simulation, which provides description of the molecular mechanisms that cause the spontaneous curvature. As a main result, for the membrane setup investigated, the spontaneous curvature is proved to be a constant plus twice the mean curvature of the membrane in its tensionless ground state. The spontaneous curvature due to the built-in transbilayer asymmetry of the membrane in terms of lipid shape is studied by the proposed method. A linear dependence of the spontaneous curvature with respect to the head-bead diameter difference and the lipid mixing ratio is discovered. The consistency with the theoretical results provides evidence supporting the validity of our method.

A coupled model, capable of simulating transonic flow, solid heat conduction, species transport, and gas radiation, is developed that provides better computational treatment of infrared radiation from hot exhaust nozzles. The modeling of gas radiation is based on a statistical narrow-band correlated-k analysis, whose parameters are deduced from the HITEMP line-by-line database. To improve computational efficiency, several methods are employed. A mixed analytical-numerical algorithm is described for the stiffness of the two-equation turbulence model and an alternating direction implicit pretreatment for the ill-conditioned matrix appearing in the coupled problem of flow and solid heat conduction. Moreover, an improved multigrid method and a symmetry plane treatment of the radiation transfer-energy equations are also introduced. Four numerical simulations are given to confirm the efficiency and accuracy of the numerical method. Finally, an account of the aerothermodynamics and infrared characteristics for two types of nozzles are presented. The infrared radiation intensity of the Chevron ejecting nozzle is clearly smaller than that of the common axisymmetric ejecting nozzle. All computations can be performed on a personal computer.

This paper presents a parallel algorithm for finding the smallest eigenvalue of a family of Hankel matrices that are ill-conditioned. Such matrices arise in random matrix theory and require the use of extremely high precision arithmetic. Surprisingly, we find that a group of commonly-used approaches that are designed for high efficiency are actually less efficient than a direct approach for this class of matrices. We then develop a parallel implementation of the algorithm that takes into account the unusually high cost of individual arithmetic operations. Our approach combines message passing and shared memory, achieving near-perfect scalability and high tolerance for network latency. We are thus able to find solutions for much larger matrices than previously possible, with the potential for extending this work to systems with greater levels of parallelism. The contributions of this work are in three areas: determination that a direct algorithm based on the secant method is more effective when extreme fixed-point precision is required than are the algorithms more typically used in parallel floating-point computations; the particular mix of optimizations required for extreme precision large matrix operations on a modern multi-core cluster, and the numerical results themselves.

This paper develops Runge-Kutta PK -based central discontinuous Galerkin (CDG) methods with WENO limiter for the one- and two-dimensional special relativistic hydrodynamical (RHD) equations, K = 1,2,3. Different from the non-central DG methods, the Runge-Kutta CDG methods have to find two approximate solutions defined on mutually dual meshes. For each mesh, the CDG approximate solutions on its dual mesh are used to calculate the flux values in the cell and on the cell boundary so that the approximate solutions on mutually dual meshes are coupled with each other, and the use of numerical flux will be avoided. The WENO limiter is adaptively implemented via two steps: the “troubled” cells are first identified by using a modified TVB minmod function, and then the WENO technique is used to locally reconstruct new polynomials of degree (2K+1) replacing the CDG solutions inside the “troubled” cells by the cell average values of the CDG solutions in the neighboring cells as well as the original cell averages of the “troubled” cells. Because the WENO limiter is only employed for finite “troubled” cells, the computational cost can be as little as possible. The accuracy of the CDG without the numerical dissipation is analyzed and calculation of the flux integrals over the cells is also addressed. Several test problems in one and two dimensions are solved by using our Runge-Kutta CDG methods with WENO limiter. The computations demonstrate that our methods are stable, accurate, and robust in solving complex RHD problems.

Recent studies on applications of the lattice Boltzmann method (LBM) and the finite-difference lattice Boltzmann method (FDLBM) to velocity slip simulations are mostly on one-dimensional (1D) problems such as a shear flow between parallel plates. Applications to a 2D problem may raise new issues. The author performed numerical simulations of rotational slip flow in coaxial cylinders as an example of 2D problem. Two types of 2D models were used. The first were multi-speed FDLBM models proposed by the author. The second was a standard LBM, the D2Q9 model. The simulations were performed applying a finite difference scheme to both the models. The study had two objectives. The first was to investigate the accuracies of LBM and FDLBM on applications to rotational slip flow. The second was to obtain an experience on application of the cylindrical coordinate system. The FDLBM model with 8 directions and the D2Q9 model showed an anisotropic flow pattern when the relaxation time constant or the Knudsen number was large. The FDLBM model with 24 directions showed accurate results even at large Knudsen numbers.

An adaptive mesh refinement strategy is proposed in this paper for the Immersed Boundary and Immersed Interface methods for two-dimensional elliptic interface problems involving singular sources. The interface is represented by the zero level set of a Lipschitz function ϕ(x,y). Our adaptive mesh refinement is done within a small tube of |ϕ(x,y)|≤δ with finer Cartesian meshes. The discrete linear system of equations is solved by a multigrid solver. The AMR methods could obtain solutions with accuracy that is similar to those on a uniform fine grid by distributing the mesh more economically, therefore, reduce the size of the linear system of the equations. Numerical examples presented show the efficiency of the grid refinement strategy.

In this paper the numerical solution of the two-dimensional sine-Gordon equation is studied. Split local artificial boundary conditions are obtained by the operator splitting method. Then the original problem is reduced to an initial boundary value problem on a bounded computational domain, which can be solved by the finite difference method. Several numerical examples are provided to demonstrate the effectiveness and accuracy of the proposed method, and some interesting propagation and collision behaviors of the solitary wave solutions are observed.

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.

The Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state is the superconducting phase for which the Cooper pairs have a non-zero total momentum, depending on the splitting of the Fermi surface sheets for electrons with opposite spin. In infinite systems the momentum is a continuous function of the temperature. In this paper, we have shown how the finite size of the system, through the discretized geometry of the Fermi surface, affects the physical properties of the FFLO state by introducing discontinuities in the Cooper pair momentum. Our calculation in an isotropic system show that the superconducting state with two opposite Cooper pair momenta is more stable than state with one momentum also in nano-size systems, where finite size effects play a crucial role.

The bubble packing method can generate high-quality node sets in simple and complex domains. However, its efficiency remains to be improved. This study is a part of an ongoing effort to introduce several acceleration schemes to reduce the cost of simulation. Firstly, allow the viscosity coefficient c in the bubble governing equations to change according the coordinate of the bubble which are defined separately as odd and normal bubbles, and meanwhile with the saw-shape relationship with time or iterations. Then, in order to relieve the over crowded initial bubble placement, two coefficients w1 and w2 are introduced to modify the insertion criterion. The range of those two coefficients are discussed to be w1 = 1, w2 ∈ [0.5,0.8]. Finally, a self-adaptive termination condition is logically set when the stable system equilibrium is achieved. Numerical examples illustrate that the computing cost can significantly decrease by roughly 80% via adopting various combination of proper schemes (except the uniform placement example), and the average qualities of corresponding Delaunay triangulation substantially exceed 0.9. It shows that those strategies are efficient and can generate a node set with high quality.

In this paper, the numerical error of two widely used methods for remapping of discrete quantities from one computational mesh to another is investigated. We compare the intuitive, but resource intensive method utilizing intersections of computational cells with the faster and simpler swept-region-based method. Both algorithms are formally second order accurate, however, they are known to produce slightly different quantity profiles in practical applications. The second-order estimate of the error formula is constructed algebraically for both algorithms so that their local accuracy can be evaluated. This general estimate is then used to assess the dependence of the performance of both methods on parameters such as the second derivatives of the remapped distribution, mesh geometry or mesh movement. Due to the complexity of such analysis, it is performed on a set of simplified elementary mesh patterns such as cell corner expansion, rotation or shear. On selected numerical tests it is demonstrated that the swept-based method can distort a symmetric quantity distribution more substantially than the intersection-based approach when the computational mesh moves in an unsuitable direction.