skip to main content


Title: A Fixed-Point Fast Sweeping WENO Method with Inverse Lax-Wendroff Boundary Treatment for Steady State of Hyperbolic Conservation Laws
Abstract

Fixed-point fast sweeping WENO methods are a class of efficient high-order numerical methods to solve steady-state solutions of hyperbolic partial differential equations (PDEs). The Gauss-Seidel iterations and alternating sweeping strategy are used to cover characteristics of hyperbolic PDEs in each sweeping order to achieve fast convergence rate to steady-state solutions. A nice property of fixed-point fast sweeping WENO methods which distinguishes them from other fast sweeping methods is that they are explicit and do not require inverse operation of nonlinear local systems. Hence, they are easy to be applied to a general hyperbolic system. To deal with the difficulties associated with numerical boundary treatment when high-order finite difference methods on a Cartesian mesh are used to solve hyperbolic PDEs on complex domains, inverse Lax-Wendroff (ILW) procedures were developed as a very effective approach in the literature. In this paper, we combine a fifth-order fixed-point fast sweeping WENO method with an ILW procedure to solve steady-state solution of hyperbolic conservation laws on complex computing regions. Numerical experiments are performed to test the method in solving various problems including the cases with the physical boundary not aligned with the grids. Numerical results show high-order accuracy and good performance of the method. Furthermore, the method is compared with the popular third-order total variation diminishing Runge-Kutta (TVD-RK3) time-marching method for steady-state computations. Numerical examples show that for most of examples, the fixed-point fast sweeping method saves more than half CPU time costs than TVD-RK3 to converge to steady-state solutions.

 
more » « less
Award ID(s):
2010107
NSF-PAR ID:
10364206
Author(s) / Creator(s):
; ; ;
Publisher / Repository:
Springer Science + Business Media
Date Published:
Journal Name:
Communications on Applied Mathematics and Computation
Volume:
5
Issue:
1
ISSN:
2096-6385
Page Range / eLocation ID:
p. 403-427
Format(s):
Medium: X
Sponsoring Org:
National Science Foundation
More Like this
  1. Embedding properties of network realizations of dissipative reduced order models Jörn Zimmerling, Mikhail Zaslavsky,Rob Remis, Shasri Moskow, Alexander Mamonov, Murthy Guddati, Vladimir Druskin, and Liliana Borcea Mathematical Sciences Department, Worcester Polytechnic Institute https://www.wpi.edu/people/vdruskin Abstract Realizations of reduced order models of passive SISO or MIMO LTI problems can be transformed to tridiagonal and block-tridiagonal forms, respectively, via dierent modications of the Lanczos algorithm. Generally, such realizations can be interpreted as ladder resistor-capacitor-inductor (RCL) networks. They gave rise to network syntheses in the rst half of the 20th century that was at the base of modern electronics design and consecutively to MOR that tremendously impacted many areas of engineering (electrical, mechanical, aerospace, etc.) by enabling ecient compression of the underlining dynamical systems. In his seminal 1950s works Krein realized that in addition to their compressing properties, network realizations can be used to embed the data back into the state space of the underlying continuum problems. In more recent works of the authors Krein's ideas gave rise to so-called nite-dierence Gaussian quadrature rules (FDGQR), allowing to approximately map the ROM state-space representation to its full order continuum counterpart on a judicially chosen grid. Thus, the state variables can be accessed directly from the transfer function without solving the full problem and even explicit knowledge of the PDE coecients in the interior, i.e., the FDGQR directly learns" the problem from its transfer function. This embedding property found applications in PDE solvers, inverse problems and unsupervised machine learning. Here we show a generalization of this approach to dissipative PDE problems, e.g., electromagnetic and acoustic wave propagation in lossy dispersive media. Potential applications include solution of inverse scattering problems in dispersive media, such as seismic exploration, radars and sonars. To x the idea, we consider a passive irreducible SISO ROM fn(s) = Xn j=1 yi s + σj , (62) assuming that all complex terms in (62) come in conjugate pairs. We will seek ladder realization of (62) as rjuj + vj − vj−1 = −shˆjuj , uj+1 − uj + ˆrj vj = −shj vj , (63) for j = 0, . . . , n with boundary conditions un+1 = 0, v1 = −1, and 4n real parameters hi, hˆi, ri and rˆi, i = 1, . . . , n, that can be considered, respectively, as the equivalent discrete inductances, capacitors and also primary and dual conductors. Alternatively, they can be viewed as respectively masses, spring stiness, primary and dual dampers of a mechanical string. Reordering variables would bring (63) into tridiagonal form, so from the spectral measure given by (62 ) the coecients of (63) can be obtained via a non-symmetric Lanczos algorithm written in J-symmetric form and fn(s) can be equivalently computed as fn(s) = u1. The cases considered in the original FDGQR correspond to either (i) real y, θ or (ii) real y and imaginary θ. Both cases are covered by the Stieltjes theorem, that yields in case (i) real positive h, hˆ and trivial r, rˆ, and in case (ii) real positive h,r and trivial hˆ,rˆ. This result allowed us a simple interpretation of (62) as the staggered nite-dierence approximation of the underlying PDE problem [2]. For PDEs in more than one variables (including topologically rich data-manifolds), a nite-dierence interpretation is obtained via a MIMO extensions in block form, e.g., [4, 3]. The main diculty of extending this approach to general passive problems is that the Stieltjes theory is no longer applicable. Moreover, the tridiagonal realization of a passive ROM transfer function (62) via the ladder network (63) cannot always be obtained in port-Hamiltonian form, i.e., the equivalent primary and dual conductors may change sign [1]. 100 Embedding of the Stieltjes problems, e.g., the case (i) was done by mapping h and hˆ into values of acoustic (or electromagnetic) impedance at grid cells, that required a special coordinate stretching (known as travel time coordinate transform) for continuous problems. Likewise, to circumvent possible non-positivity of conductors for the non-Stieltjes case, we introduce an additional complex s-dependent coordinate stretching, vanishing as s → ∞ [1]. This stretching applied in the discrete setting induces a diagonal factorization, removes oscillating coecients, and leads to an accurate embedding for moderate variations of the coecients of the continuum problems, i.e., it maps discrete coecients onto the values of their continuum counterparts. Not only does this embedding yields an approximate linear algebraic algorithm for the solution of the inverse problems for dissipative PDEs, it also leads to new insight into the properties of their ROM realizations. We will also discuss another approach to embedding, based on Krein-Nudelman theory [5], that results in special data-driven adaptive grids. References [1] Borcea, Liliana and Druskin, Vladimir and Zimmerling, Jörn, A reduced order model approach to inverse scattering in lossy layered media, Journal of Scientic Computing, V. 89, N1, pp. 136,2021 [2] Druskin, Vladimir and Knizhnerman, Leonid, Gaussian spectral rules for the three-point second dierences: I. A two-point positive denite problem in a semi-innite domain, SIAM Journal on Numerical Analysis, V. 37, N 2, pp.403422, 1999 [3] Druskin, Vladimir and Mamonov, Alexander V and Zaslavsky, Mikhail, Distance preserving model order reduction of graph-Laplacians and cluster analysis, Druskin, Vladimir and Mamonov, Alexander V and Zaslavsky, Mikhail, Journal of Scientic Computing, V. 90, N 1, pp 130, 2022 [4] Druskin, Vladimir and Moskow, Shari and Zaslavsky, Mikhail LippmannSchwingerLanczos algorithm for inverse scattering problems, Inverse Problems, V. 37, N. 7, 2021, [5] Mark Adolfovich Nudelman The Krein String and Characteristic Functions of Maximal Dissipative Operators, Journal of Mathematical Sciences, 2004, V 124, pp 49184934 Go back to Plenary Speakers Go back to Speakers Go back 
    more » « less
  2. Abstract In this article, the recently discovered phenomenon of delayed Hopf bifurcations (DHB) in reaction–diffusion partial differential equations (PDEs) is analysed in the cubic Complex Ginzburg–Landau equation, as an equation in its own right, with a slowly varying parameter. We begin by using the classical asymptotic methods of stationary phase and steepest descents on the linearized PDE to show that solutions, which have approached the attracting quasi-steady state (QSS) before the Hopf bifurcation remain near that state for long times after the instantaneous Hopf bifurcation and the QSS has become repelling. In the complex time plane, the phase function of the linearized PDE has a saddle point, and the Stokes and anti-Stokes lines are central to the asymptotics. The non-linear terms are treated by applying an iterative method to the mild form of the PDE given by perturbations about the linear particular solution. This tracks the closeness of solutions near the attracting and repelling QSS in the full, non-linear PDE. Next, we show that beyond a key Stokes line through the saddle there is a curve in the space-time plane along which the particular solution of the linear PDE ceases to be exponentially small, causing the solution of the non-linear PDE to diverge from the repelling QSS and exhibit large-amplitude oscillations. This curve is called the space–time buffer curve. The homogeneous solution also stops being exponentially small in a spatially dependent manner, as determined also by the initial data and time. Hence, a competition arises between these two solutions, as to which one ceases to be exponentially small first, and this competition governs spatial dependence of the DHB. We find four different cases of DHB, depending on the outcomes of the competition, and we quantify to leading order how these depend on the main system parameters, including the Hopf frequency, initial time, initial data, source terms, and diffusivity. Examples are presented for each case, with source terms that are a uni-modal function, a smooth step function, a spatially periodic function and an algebraically growing function. Also, rich spatio-temporal dynamics are observed in the post-DHB oscillations. Finally, it is shown that large-amplitude source terms can be designed so that solutions spend substantially longer times near the repelling QSS, and hence, region-specific control over the delayed onset of oscillations can be achieved. 
    more » « less
  3. Chi-Wang Shu (Ed.)
    GPU computing is expected to play an integral part in all modern Exascale supercomputers. It is also expected that higher order Godunov schemes will make up about a significant fraction of the application mix on such supercomputers. It is, therefore, very important to prepare the community of users of higher order schemes for hyperbolic PDEs for this emerging opportunity. Not every algorithm that is used in the space-time update of the solution of hyperbolic PDEs will take well to GPUs. However, we identify a small core of algorithms that take exceptionally well to GPU computing. Based on an analysis of available options, we have been able to identify weighted essentially non-oscillatory (WENO) algorithms for spatial reconstruction along with arbitrary derivative (ADER) algorithms for time extension followed by a corrector step as the winning three-part algorithmic combination. Even when a winning subset of algorithms has been identified, it is not clear that they will port seamlessly to GPUs. The low data throughput between CPU and GPU, as well as the very small cache sizes on modern GPUs, implies that we have to think through all aspects of the task of porting an application to GPUs. For that reason, this paper identifies the techniques and tricks needed for making a successful port of this very useful class of higher order algorithms to GPUs. Application codes face a further challenge—the GPU results need to be practically indistinguishable from the CPU results—in order for the legacy knowledge bases embedded in these applications codes to be preserved during the port of GPUs. This requirement often makes a complete code rewrite impossible. For that reason, it is safest to use an approach based on OpenACC directives, so that most of the code remains intact (as long as it was originally well-written). This paper is intended to be a one-stop shop for anyone seeking to make an OpenACC-based port of a higher order Godunov scheme to GPUs. We focus on three broad and high-impact areas where higher order Godunov schemes are used. The first area is computational fluid dynamics (CFD). The second is computational magnetohydrodynamics (MHD) which has an involution constraint that has to be mimetically preserved. The third is computational electrodynamics (CED) which has involution constraints and also extremely stiff source terms. Together, these three diverse uses of higher order Godunov methodology, cover many of the most important applications areas. In all three cases, we show that the optimal use of algorithms, techniques, and tricks, along with the use of OpenACC, yields superlative speedups on GPUs. As a bonus, we find a most remarkable and desirable result: some higher order schemes, with their larger operations count per zone, show better speedup than lower order schemes on GPUs. In other words, the GPU is an optimal stratagem for overcoming the higher computational complexities of higher order schemes. Several avenues for future improvement have also been identified. A scalability study is presented for a real-world application using GPUs and comparable numbers of high-end multicore CPUs. It is found that GPUs offer a substantial performance benefit over comparable number of CPUs, especially when all the methods designed in this paper are used. 
    more » « less
  4. Abstract

    Solving analytically intractable partial differential equations (PDEs) that involve at least one variable defined on an unbounded domain arises in numerous physical applications. Accurately solving unbounded domain PDEs requires efficient numerical methods that can resolve the dependence of the PDE on the unbounded variable over at least several orders of magnitude. We propose a solution to such problems by combining two classes of numerical methods: (i) adaptive spectral methods and (ii) physics-informed neural networks (PINNs). The numerical approach that we develop takes advantage of the ability of PINNs to easily implement high-order numerical schemes to efficiently solve PDEs and extrapolate numerical solutions at any point in space and time. We then show how recently introduced adaptive techniques for spectral methods can be integrated into PINN-based PDE solvers to obtain numerical solutions of unbounded domain problems that cannot be efficiently approximated by standard PINNs. Through a number of examples, we demonstrate the advantages of the proposed spectrally adapted PINNs in solving PDEs and estimating model parameters from noisy observations in unbounded domains.

     
    more » « less
  5. null (Ed.)
    Abstract We outline and interpret a recently developed theory of impedance matching or reflectionless excitation of arbitrary finite photonic structures in any dimension. The theory includes both the case of guided wave and free-space excitation. It describes the necessary and sufficient conditions for perfectly reflectionless excitation to be possible and specifies how many physical parameters must be tuned to achieve this. In the absence of geometric symmetries, such as parity and time-reversal, the product of parity and time-reversal, or rotational symmetry, the tuning of at least one structural parameter will be necessary to achieve reflectionless excitation. The theory employs a recently identified set of complex frequency solutions of the Maxwell equations as a starting point, which are defined by having zero reflection into a chosen set of input channels, and which are referred to as R-zeros. Tuning is generically necessary in order to move an R-zero to the real frequency axis, where it becomes a physical steady-state impedance-matched solution, which we refer to as a reflectionless scattering mode (RSM). In addition, except in single-channel systems, the RSM corresponds to a particular input wavefront, and any other wavefront will generally not be reflectionless. It is useful to consider the theory as representing a generalization of the concept of critical coupling of a resonator, but it holds in arbitrary dimension, for arbitrary number of channels, and even when resonances are not spectrally isolated. In a structure with parity and time-reversal symmetry (a real dielectric function) or with parity–time symmetry, generically a subset of the R-zeros has real frequencies, and reflectionless states exist at discrete frequencies without tuning. However, they do not exist within every spectral range, as they do in the special case of the Fabry–Pérot or two-mirror resonator, due to a spontaneous symmetry-breaking phenomenon when two RSMs meet. Such symmetry-breaking transitions correspond to a new kind of exceptional point, only recently identified, at which the shape of the reflection and transmission resonance lineshape is flattened. Numerical examples of RSMs are given for one-dimensional multimirror cavities, a two-dimensional multiwaveguide junction, and a multimode waveguide functioning as a perfect mode converter. Two solution methods to find R-zeros and RSMs are discussed. The first one is a straightforward generalization of the complex scaling or perfectly matched layer method and is applicable in a number of important cases; the second one involves a mode-specific boundary matching method that has only recently been demonstrated and can be applied to all geometries for which the theory is valid, including free space and multimode waveguide problems of the type solved here. 
    more » « less