Journal of Computational Physics: X最新文献

In this paper, we develop sparse grid discontinuous Galerkin (DG) schemes for the Vlasov-Maxwell (VM) equations. The VM system is a fundamental kinetic model in plasma physics, and its numerical computations are quite demanding, due to its intrinsic high-dimensionality and the need to retain many properties of the physical solutions. To break the curse of dimensionality, we consider the sparse grid DG methods that were recently developed in [20], [21] for transport equations. Such methods are based on multiwavelets on tensorized nested grids and can significantly reduce the numbers of degrees of freedom. We formulate two versions of the schemes: sparse grid DG and adaptive sparse grid DG methods for the VM system. Their key properties and implementation details are discussed. Accuracy and robustness are demonstrated by numerical tests, with emphasis on comparison of the performance of the two methods, as well as with their full grid counterparts.

In this paper we present and analyze finite difference numerical schemes for the Cahn-Hilliard equation with a logarithmic Flory Huggins energy potential. Both first and second order accurate temporal algorithms are considered. In the first order scheme, we treat the nonlinear logarithmic terms and the surface diffusion term implicitly, and update the linear expansive term and the mobility explicitly. We provide a theoretical justification that this numerical algorithm has a unique solution, such that the positivity is always preserved for the logarithmic arguments, i.e., the phase variable is always between −1 and 1, at a point-wise level. In particular, our analysis reveals a subtle fact: the singular nature of the logarithmic term around the values of −1 and 1 prevents the numerical solution reaching these singular values, so that the numerical scheme is always well-defined as long as the numerical solution stays similarly bounded at the previous time step. Furthermore, an unconditional energy stability of the numerical scheme is derived, without any restriction for the time step size. Such an analysis technique can also be applied to a second order numerical scheme in which the BDF temporal stencil is applied, the expansive term is updated by a second order Adams-Bashforth explicit extrapolation formula, and an artificial Douglas-Dupont regularization term is added to ensure the energy dissipativity. The unique solvability and the positivity-preserving property for the second order scheme are proved using similar ideas, namely, the singular nature of the logarithmic term plays an essential role. For both the first and second order accurate schemes, we are able to derive an optimal rate convergence analysis. The case with a non-constant mobility is analyzed as well. We also describe a practical and efficient multigrid solver for the proposed numerical schemes, and present some numerical results, which demonstrate the robustness of the numerical schemes.

We implemented a fast Reciprocal Monte Carlo algorithm to accurately solve radiative heat transfer in turbulent flows of non-grey participating media that can be coupled to fully resolved turbulent flows, namely to Direct Numerical Simulation (DNS). The spectrally varying absorption coefficient is treated in a narrow-band fashion with a correlated-k distribution. The implementation is verified with analytical solutions and validated with results from literature and line-by-line Monte Carlo computations. The method is implemented on GPU with a thorough attention to memory transfer and computational efficiency. The bottlenecks that dominate the computational expenses are addressed, and several techniques are proposed to optimize the GPU execution. By implementing the proposed algorithmic accelerations, while maintaining the same accuracy, a speed-up of up to 3 orders of magnitude can be achieved.

The Euler equations for inviscid and compressible flows are used for modeling interfacial instabilities, and in doing so all the physical dissipations are ignored under the consideration that they are extremely weak. However, numerical simulations of interfacial instabilities with numerical dissipations or with little dissipations suffer from nonphysical artifacts on the interfaces in late times of the interfacial developments. In this paper we introduce numerical dissipations for our previously developed conservative front-tracking method that simulate tangentially the missing physical dissipations in the Euler equations on the interfaces. Numerical examples show that they suppress numerical artifacts on the tracked interfaces and help to accomplish the simulations of interfacial instabilities on fine grids.

We present the results of large scale simulations of 4th order nonlinear partial differential equations of diffusion type that are typically encountered when modeling dynamics of thin fluid films on substrates. The simulations are based on the alternate direction implicit (ADI) method, with the main part of the computational work carried out in the GPU computing environment. Efficient and accurate computations allow for simulations on large computational domains in three spatial dimensions (3D) and for long computational times. We apply the methods developed to the particular problem of instabilities of thin fluid films of nanoscale thickness. The large scale of the simulations minimizes the effects of boundaries, and also allows for simulating domains of the size encountered in published experiments. As an outcome, we can analyze the development of instabilities with an unprecedented level of detail. A particular focus is on analyzing the manner in which instability develops, in particular regarding differences between spinodal and nucleation types of dewetting for linearly unstable films, as well as instabilities of metastable films. Simulations in 3D allow for consideration of some recent results that were previously obtained in the 2D geometry [28]. Some of the new results include using Fourier transforms as well as topological invariants (Betti numbers) to distinguish the outcomes of spinodal and nucleation types of instabilities, describing in precise terms the complex processes that lead to the formation of satellite drops, as well as distinguishing the shape of the evolving film front in linearly unstable and metastable regimes. We also discuss direct comparison between simulations and available experimental results for nematic liquid crystal and polymer films.

A significant drawback of Lagrangian (particle-tracking) reactive transport models has been their inability to properly simulate interactions between solid and liquid chemical phases, such as dissolution and precipitation reactions. This work addresses that problem by implementing a mass-transfer algorithm between mobile and immobile sets of particles that allows aqueous species of reactant that are undergoing transport to interact with stationary solid species. This mass-transfer algorithm is demonstrated to solve the diffusion equation for an arbitrarily small level of diffusion and thus does not introduce any spurious mixing. The algorithm can be combined with random walks to simulate the desired total level of diffusion in a reactive transport system.

The Du Fort–Frankel scheme for the one-dimensional Schrödinger equation is shown to be equivalent, under a time-dependent unitary transformation, to the Ablowitz–Kruskal–Ladik scheme for the Klein–Gordon equation. The Schrödinger equation describes a non-relativistic quantum particle, while the Klein–Gordon equation describes a relativistic particle. The conditional convergence of the Du Fort–Frankel scheme to solutions of the Schrödinger equation arises because solutions of the Klein–Gordon equation only approximate solutions of the Schrödinger equation in the non-relativistic limit. The time-dependent unitary transformation is the discrete analog of the transformation that arises from seeking a non-relativistic limit using the interaction picture of quantum mechanics to decompose the Klein–Gordon Hamiltonian into the relativistic rest energy and a remainder. The Ablowitz–Kruskal–Ladik scheme is in turn decomposed into a quantum lattice gas automaton for the one-dimensional Dirac equation, which is also the one-dimensional discrete time quantum walk. This relativistic interpretation clarifies the origin of the known discrete invariant of the Du Fort–Frankel scheme as expressing conservation of probability for the 2-component wavefunction in the one-dimensional Dirac equation under discrete unitary evolution. It also leads to a second invariant, the matrix element of the evolution operator, whose imaginary part gives a discrete approximation to the expectation of the non-relativistic Schrödinger Hamiltonian.

We propose a new numerical method to solve the Cahn-Hilliard equation coupled with non-linear wetting boundary conditions. We show that the method is mass-conservative and that the discrete solution satisfies a discrete energy law similar to the one satisfied by the exact solution. We perform several tests inspired by realistic situations to verify the accuracy and performance of the method: wetting of a chemically heterogeneous substrate in three dimensions, wetting-driven nucleation in a complex two-dimensional domain and three-dimensional diffusion through a porous medium.

We present a new method for evolving the equations of magnetohydrodynamics (both Newtonian and relativistic) that is capable of maintaining a divergence-free magnetic field ($\nabla \cdot B=0$) on adaptively refined, conformally moving meshes. The method relies on evolving the magnetic vector potential and then using it to reconstruct the magnetic fields. The advantage of this approach is that the vector potential is not subject to a constraint equation in the same way the magnetic field is, and so can be refined and moved in a straightforward way. We test this new method against a wide array of problems from simple Alfvén waves on a uniform grid to general relativistic MHD simulations of black hole accretion on a nested, spherical-polar grid. We find that the code produces accurate results and in all cases maintains a divergence-free magnetic field to machine precision.

We present a parallel hp-adaptive high order (spectral) discontinuous Galerkin method for approximation of the incompressible Navier-Stokes equations. The spatial discretization consists of equal-order polynomial approximations of the fluid velocity and pressure via discontinuous Galerkin spatial discretizations. For the nonlinear convective term we select the local Lax-Friedrichs flux, while for the divergence and gradient operators central fluxes are chosen. For the diffusive term, we use an interior penalty discontinuous Galerkin method to ensure stability and invertibility. The temporal discretization is an implicit-explicit Runge-Kutta method paired with a high-order splitting procedure to efficiently enforce the incompressibility condition at each time step. The compact stencil size, explicit time stepping of nonlinear terms, and inversion of sparse linear systems make the resulting method simple to parallelize while the local nature of the discontinuous Galerkin approximation makes hp-adaptive refinement natural to implement. We detail our implementation consisting of a tensor product basis of high order polynomials on quadrilateral elements, and implement hp-adaptivity using an inexpensive a posteriori error estimator to determine where refinement is necessary. p-Multigrid and pressure projection techniques are used to precondition the conjugate gradient linear solvers. We present several numerical tests to demonstrate the efficacy of the method, in particular in reducing the number of degrees of freedom needed and allocating computing resources to regions of sharp variation in transient incompressible Navier-Stokes flows.