skip to main content


Title: A mass- and energy-conserved DG method for the Schrödinger-Poisson equation
In this paper, we construct, analyze, and numerically validate conservative discontinuous Galerkin (DG) schemes for approximating the Schr\"{o}dinger-Poisson equation. The proposed schemes all satisfy both mass and energy conservation. For the semi-discrete DG scheme optimal $L^2$ error estimates are obtained. Efficient iterative solvers are also constructed to solve the second order implicit time discretization. A number of numerical tests are presented to demonstrate the method’s accuracy and robustness, confirming that both mass and energy are well preserved over long time simulations.  more » « less
Award ID(s):
1812666
NSF-PAR ID:
10282820
Author(s) / Creator(s):
;
Date Published:
Journal Name:
Numerical Algorithms
ISSN:
1017-1398
Format(s):
Medium: X
Sponsoring Org:
National Science Foundation
More Like this
  1. In this paper, we consider Maxwell’s equations in linear dispersive media described by a single-pole Lorentz model for electronic polarization. We study two classes of commonly used spatial discretizations: finite difference methods (FD) with arbitrary even order accuracy in space and high spatial order discontinuous Galerkin (DG) finite element methods. Both types of spatial discretizations are coupled with second order semi-implicit leap-frog and implicit trapezoidal temporal schemes. By performing detailed dispersion analysis for the semi-discrete and fully discrete schemes, we obtain rigorous quantification of the dispersion error for Lorentz dispersive dielectrics. In particular, comparisons of dispersion error can be made taking into account the model parameters, and mesh sizes in the design of the two types of schemes. This work is a continuation of our previous research on energy-stable numerical schemes for nonlinear dispersive optical media [6,7]. The results for the numerical dispersion analysis of the reduced linear model, considered in the present paper, can guide us in the optimal choice of discretization parameters for the more complicated and nonlinear models. The numerical dispersion analysis of the fully discrete FD and DG schemes, for the dispersive Maxwell model considered in this paper, clearly indicate the dependence of the numerical dispersion errors on spatial and temporal discretizations, their order of accuracy, mesh discretization parameters and model parameters. The results obtained here cannot be arrived at by considering discretizations of Maxwell’s equations in free space. In particular, our results contrast the advantages and disadvantages of using high order FD or DG schemes and leap-frog or trapezoidal time integrators over different frequency ranges using a variety of measures 
    more » « less
  2. Abstract

    This paper examines a class of involution-constrained PDEs where some part of the PDE system evolves a vector field whose curl remains zero or grows in proportion to specified source terms. Such PDEs are referred to as curl-free or curl-preserving, respectively. They arise very frequently in equations for hyperelasticity and compressible multiphase flow, in certain formulations of general relativity and in the numerical solution of Schrödinger’s equation. Experience has shown that if nothing special is done to account for the curl-preserving vector field, it can blow up in a finite amount of simulation time. In this paper, we catalogue a class of DG-like schemes for such PDEs. To retain the globally curl-free or curl-preserving constraints, the components of the vector field, as well as their higher moments, must be collocated at the edges of the mesh. They are updated using potentials collocated at the vertices of the mesh. The resulting schemes: (i) do not blow up even after very long integration times, (ii) do not need any special cleaning treatment, (iii) can operate with large explicit timesteps, (iv) do not require the solution of an elliptic system and (v) can be extended to higher orders using DG-like methods. The methods rely on a special curl-preserving reconstruction and they also rely on multidimensional upwinding. The Galerkin projection, highly crucial to the design of a DG method, is now conducted at the edges of the mesh and yields a weak form update that uses potentials obtained at the vertices of the mesh with the help of a multidimensional Riemann solver. A von Neumann stability analysis of the curl-preserving methods is conducted and the limiting CFL numbers of this entire family of methods are catalogued in this work. The stability analysis confirms that with the increasing order of accuracy, our novel curl-free methods have superlative phase accuracy while substantially reducing dissipation. We also show that PNPM-like methods, which only evolve the lower moments while reconstructing the higher moments, retain much of the excellent wave propagation characteristics of the DG-like methods while offering a much larger CFL number and lower computational complexity. The quadratic energy preservation of these methods is also shown to be excellent, especially at higher orders. The methods are also shown to be curl-preserving over long integration times.

     
    more » « less
  3. Abstract

    This article presents high order accurate discontinuous Galerkin (DG) methods for wave problems on moving curved meshes with general choices of basis and quadrature. The proposed method adopts an arbitrary Lagrangian–Eulerian formulation to map the wave equation from a time‐dependent moving physical domain onto a fixed reference domain. For moving curved meshes, weighted mass matrices must be assembled and inverted at each time step when using explicit time‐stepping methods. We avoid this step by utilizing an easily invertible weight‐adjusted approximation. The resulting semi‐discrete weight‐adjusted DG scheme is provably energy stable up to a term that (for a fixed time interval) converges to zero with the same rate as the optimal error estimate. Numerical experiments using both polynomial and B‐spline bases verify the high order accuracy and energy stability of proposed methods.

     
    more » « less
  4. In this article, we introduce an invariant‐region‐preserving (IRP) limiter for thep‐system and the corresponding viscousp‐system, both of which share the same invariant region. Rigorous analysis is presented to show that for smooth solutions the order of approximation accuracy is not destroyed by the IRP limiter, provided the cell average stays away from the boundary of invariant region. Moreover, this limiter is explicit, and easy for computer implementation. A generic algorithm incorporating the IRP limiter is presented for high order finite volume type schemes as long as the evolved cell average of the underlying scheme stays strictly within the invariant region. For high order discontinuous Galerkin (DG) schemes to thep‐system, sufficient conditions are obtained for cell averages to stay in the invariant region. For the viscousp‐system, we design both second and third order IRP DG schemes. Numerical experiments are provided to test the proven properties of the IRP limiter and the performance of IRP DG schemes.

     
    more » « less
  5. Sub-cloud rain evaporation in the trade wind region significantly influences the boundary layer mass and energy budgets. Parameterizing it is, however, difficult due to the sparsity of well-resolved rain observations and the challenges of sampling short-lived marine cumulus clouds. In this study, sub-cloud rain evaporation is analyzed using a steady-state, one-dimensional model that simulates changes in drop sizes, relative humidity, and rain isotopic composition. The model is initialized with relative humidity, raindrop size distributions, and water vapor isotope ratios (e.g., δDv, δ18Ov) sampled by the NOAA P3 aircraft during the Atlantic Tradewind Ocean–Atmosphere Mesoscale Interaction Campaign (ATOMIC), which was part of the larger EUREC4A (ElUcidating the RolE of Clouds–Circulation Coupling in ClimAte) field program. The modeled surface precipitation isotope ratios closely match the observations from EUREC4A ground-based and ship-based platforms, lending credibility to our model. The model suggests that 63 % of the rain mass evaporates in the sub-cloud layer across 22 P3 cases. The vertical distribution of the evaporated rain flux is top heavy for a narrow (σ) raindrop size distribution (RSD) centered over a small geometric mean diameter (Dg) at the cloud base. A top-heavy profile has a higher rain-evaporated fraction (REF) and larger changes in the rain deuterium excess (d=δD-8×δ18O) between the cloud base and the surface than a bottom-heavy profile, which results from a wider RSD with larger Dg. The modeled REF and change in d are also more strongly influenced by cloud base Dg and σ rather than the concentration of raindrops. The model results are accurate as long as the variations in the relative humidity conditions are accounted for. Relative humidity alone, however, is a poor indicator of sub-cloud rain evaporation. Overall, our analysis indicates the intricate dependence of sub-cloud rain evaporation on both thermodynamic and microphysical processes in the trade wind region. 
    more » « less