skip to main content

Title: High Spatial Order Energy Stable FDTD Methods for Maxwell’s Equations in Nonlinear Optical Media in One Dimension
In this paper, we consider electromagnetic (EM) wave propagation in nonlinear optical media in one spatial dimension. We model the EM wave propagation by the time- dependent Maxwell’s equations coupled with a system of nonlinear ordinary differential equations (ODEs) for the response of the medium to the EM waves. The nonlinearity in the ODEs describes the instantaneous electronic Kerr response and the residual Raman molecular vibrational response. The ODEs also include the single resonance linear Lorentz dispersion. For such model, we will design and analyze fully discrete finite difference time domain (FDTD) methods that have arbitrary (even) order in space and second order in time. It is challenging to achieve provable stability for fully discrete methods, and this depends on the choices of temporal discretizations of the nonlinear terms. In Bokil et al. (J Comput Phys 350:420–452, 2017), we proposed novel modifications of second-order leap-frog and trapezoidal temporal schemes in the context of discontinuous Galerkin methods to discretize the nonlinear terms in this Maxwell model. Here, we continue this work by developing similar time discretizations within the framework of FDTD methods. More specifically, we design fully discrete modified leap-frog FDTD methods which are proved to be stable under appropriate CFL conditions. These method can be viewed as an extension of the Yee-FDTD scheme to this more » nonlinear Maxwell model. We also design fully discrete trapezoidal FDTD methods which are proved to be unconditionally stable. The performance of the fully discrete FDTD methods are demonstrated through numerical experiments involving kink, antikink waves and third harmonic generation in soliton propagation. « less
; ; ;
Award ID(s):
Publication Date:
Journal Name:
Journal of scientific computing
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 dispersionmore »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« less
  2. In this paper, we apply two fully-discrete local discontinuous Galerkin (LDG) methods to the compressible wormhole propagation. We will prove the stability and error estimates of the schemes. Traditional LDG methods use the diffusion term to control of convection term to obtain the stability for some linear equations. However, the variables in wormhole propagation are coupled together and the whole system is highly nonlinear. Therefore, it is extremely difficult to obtain the stability for fully-discrete LDG methods. To fix this gap, we introduce a new auxiliary variable including both the convection and diffusion terms. Moreover, we also construct a special time integration for the porosity, leading to physically relevant numerical approximations and controllable growth rate of the porosity. With a reasonable growth rate, it is possible to handle the time level mismatch in the first-order fully discrete scheme and obtain the stability of the scheme. For the whole system, we will prove that under weak temporal-spatial conditions, the optimal error estimates for the pressure, velocity, porosity and concentration under different norms can be obtained. Numerical experiments are also given to verify the theoretical results.
  3. The thermal radiative transfer (TRT) equations form an integro-differential system that describes the propagation and collisional interactions of photons. Computing accurate and efficient numerical solutions TRT are challenging for several reasons, the first of which is that TRT is defined on a high-dimensional phase space that includes the independent variables of time, space, and velocity. In order to reduce the dimensionality of the phase space, classical approaches such as the P$_N$ (spherical harmonics) or the S$_N$ (discrete ordinates) ansatz are often used in the literature. In this work, we introduce a novel approach: the hybrid discrete (H$^T_N$) approximation to the radiative thermal transfer equations. This approach acquires desirable properties of both P$_N$ and S$_N$, and indeed reduces to each of these approximations in various limits: H$^1_N$ $\equiv$ P$_N$ and H$^T_0$ $\equiv$ S$_T$. We prove that H$^T_N$ results in a system of hyperbolic partial differential equations for all $T\ge 1$ and $N\ge 0$. Another challenge in solving the TRT system is the inherent stiffness due to the large timescale separation between propagation and collisions, especially in the diffusive (i.e., highly collisional) regime. This stiffness challenge can be partially overcome via implicit time integration, although fully implicit methods may become computationally expensivemore »due to the strong nonlinearity and system size. On the other hand, explicit time-stepping schemes that are not also asymptotic-preserving in the highly collisional limit require resolving the mean-free path between collisions, making such schemes prohibitively expensive. In this work we develop a numerical method that is based on a nodal discontinuous Galerkin discretization in space, coupled with a semi-implicit discretization in time. In particular, we make use of a second order explicit Runge-Kutta scheme for the streaming term and an implicit Euler scheme for the material coupling term. Furthermore, in order to solve the material energy equation implicitly after each predictor and corrector step, we linearize the temperature term using a Taylor expansion; this avoids the need for an iterative procedure, and therefore improves efficiency. In order to reduce unphysical oscillation, we apply a slope limiter after each time step. Finally, we conduct several numerical experiments to verify the accuracy, efficiency, and robustness of the H$^T_N$ ansatz and the numerical discretizations.« less

    Tsunami generation by offshore earthquakes is a problem of scientific interest and practical relevance, and one that requires numerical modelling for data interpretation and hazard assessment. Most numerical models utilize two-step methods with one-way coupling between separate earthquake and tsunami models, based on approximations that might limit the applicability and accuracy of the resulting solution. In particular, standard methods focus exclusively on tsunami wave modelling, neglecting larger amplitude ocean acoustic and seismic waves that are superimposed on tsunami waves in the source region. In this study, we compare four earthquake-tsunami modelling methods. We identify dimensionless parameters to quantitatively approximate dominant wave modes in the earthquake-tsunami source region, highlighting how the method assumptions affect the results and discuss which methods are appropriate for various applications such as interpretation of data from offshore instruments in the source region. Most methods couple a 3-D solid earth model, which provides the seismic wavefield or at least the static elastic displacements, with a 2-D depth-averaged shallow water tsunami model. Assuming the ocean is incompressible and tsunami propagation is negligible over the earthquake duration leads to the instantaneous source method, which equates the static earthquake seafloor uplift with the initial tsunami sea surface height. Formore »longer duration earthquakes, it is appropriate to follow the time-dependent source method, which uses time-dependent earthquake seafloor velocity as a forcing term in the tsunami mass balance. Neither method captures ocean acoustic or seismic waves, motivating more advanced methods that capture the full wavefield. The superposition method of Saito et al. solves the 3-D elastic and acoustic equations to model the seismic wavefield and response of a compressible ocean without gravity. Then, changes in sea surface height from the zero-gravity solution are used as a forcing term in a separate tsunami simulation, typically run with a shallow water solver. A superposition of the earthquake and tsunami solutions provides an approximation to the complete wavefield. This method is algorithmically a two-step method. The complete wavefield is captured in the fully coupled method, which utilizes a coupled solid Earth and compressible ocean model with gravity. The fully coupled method, recently incorporated into the 3-D open-source code SeisSol, simultaneously solves earthquake rupture, seismic waves and ocean response (including gravity). We show that the superposition method emerges as an approximation to the fully coupled method subject to often well-justified assumptions. Furthermore, using the fully coupled method, we examine how the source spectrum and ocean depth influence the expression of oceanic Rayleigh waves. Understanding the range of validity of each method, as well as its computational expense, facilitates the selection of modelling methods for the accurate assessment of earthquake and tsunami hazards and the interpretation of data from offshore instruments.

    « less
  5. Abstract Adaptive time stepping methods for metastable dynamics of the Allen–Cahn and Cahn–Hilliard equations are investigated in the spatially continuous, semi-discrete setting. We analyse the performance of a number of first and second order methods, formally predicting step sizes required to satisfy specified local truncation error $$\sigma $$ σ in the limit of small length scale parameter $$\epsilon \rightarrow 0$$ ϵ → 0 during meta-stable dynamics. The formal predictions are made under stability assumptions that include the preservation of the asymptotic structure of the diffuse interface, a concept we call profile fidelity. In this setting, definite statements about the relative behaviour of time stepping methods can be made. Some methods, including all so-called energy stable methods but also some fully implicit methods, require asymptotically more time steps than others. The formal analysis is confirmed in computational studies. We observe that some provably energy stable methods popular in the literature perform worse than some more standard schemes. We show further that when Backward Euler is applied to meta-stable Allen–Cahn dynamics, the energy decay and profile fidelity properties for these discretizations are preserved for much larger time steps than previous analysis would suggest. The results are established asymptotically for general interfaces, withmore »a rigorous proof for radial interfaces. It is shown analytically and computationally that for most reaction terms, Eyre type time stepping performs asymptotically worse due to loss of profile fidelity.« less