Local Time Stepping Research Articles

A modified dual time integration technique for aerodynamic quasi-static and dynamic stall hysteresis

Simulation of the aerodynamic stall phenomenon in both quasi-static and dynamic conditions requires expensive computational resources. The computations become even more costly for static stall hysteresis using an unsteady solver with very slow variation of angle of attack at low reduced frequencies. In an explicit time-marching solver that satisfies the low Courant number condition, that is, [Formula: see text], the computational cost for such simulations becomes prohibitive, especially at higher Reynolds numbers due to the presence of thin-stretched cells with large aspect ratio in the boundary layer. In this paper, a segregated solver method such as the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) is modified as a dual pseudo-time marching method so that the unsteady problem at each time step is reformulated as a steady state problem. The resulting system of equations in the discretized finite volume formulation is then reduced to zero or near-zero residuals using available convergence acceleration methods such as local time stepping, multi-grid acceleration and residual smoothing. The performance and accuracy of the implemented algorithm was tested for three different airfoils at low to moderate Reynolds numbers in both incompressible and compressible flow conditions covering both attached and separated flow regimes. The results obtained are in close agreement with the published experimental and computational results for both quasi-static and dynamic stall and have demonstrated significant savings in computational time.

A Minimal Round-Trip Strategy Based on Graph Matching for Parallel DGTD Method With Local Time-Stepping

The load distribution in the local time stepping (LTS) method significantly impacts its computing efficiency. This letter proposes a minimal round-trip (MRT) strategy of the LTS method to balance the communication load of the discontinuous Galerkin time-domain (DGTD) method. By discovering the matching of the connected graph of computing nodes, independent communication with a similar load can be done in the same round trip to minimize waiting between nodes in nonblocking communication, thereby decreasing the communication time of the DGTD-LTS technique. The numerical results indicate that the MRT strategy reduces the communication time between processors by 50% and improves the parallel performance when the LTS method is implemented in the DGTD method. The parallel scale of the MRT approach may be increased to 16 000 nodes (1 040 000 cores) on the supercomputer, and the parallel efficiency is greater than 73.8%.

Storm Surge Modeling as an Application of Local Time‐Stepping in MPAS‐Ocean

AbstractThis paper presents the first practical application of local time‐stepping (LTS) schemes in the Model for Prediction Across Scales‐Ocean (MPAS‐O). We use LTS schemes in a single‐layer, global ocean model that predicts the storm surge around the eastern coast of the United States during Hurricane Sandy. The variable‐resolution meshes used are of unprecedentedly high resolution in MPAS‐O, containing cells as small as 125 m wide in Delaware Bay. It is shown that a particular, third‐order LTS scheme (LTS3) produces sea‐surface height solutions that are of comparable quality to solutions produced by the classical four‐stage, fourth‐order Runge‐Kutta method (RK4) with a uniform time step on the same meshes. Furthermore, LTS3 is up to 35% faster in the best cases considered, where the number of cells using the coarse time‐step relative to those using the fine time‐step is as low as 1:1. This shows that LTS schemes are viable for use in MPAS‐O with the added benefit of substantially less computational cost. The results of these performance experiments inform us of the requirements for efficient mesh design and configuration of LTS regions for LTS schemes. In particular, we see that for LTS to be efficient on a given mesh, it is important to have enough cells using the coarse time‐step relative to those using the fine time‐step, typically at least 1:5 to see an increase in performance.

A robust discontinuous Galerkin time domain algorithm for penetrable thin layers

<sec>The thin layer problem is common in some engineering cases, such as a coating target or dielectric radome. If the skin depth of the medium behind the thin layer is much smaller than the medium scale, the transmission line impedance boundary can be used to solve the thin layer problem . When both sides of the thin layer are free space or the skin depth of the medium behind the thin layer is greater than or comparable to the scale of the medium, the transmission line impedance boundary condition is not applicable. Such problems need incident and reflected wave data and are also known as penetrable thin-layer problems.</sec><sec>When dealing with penetrable thin-layer problems, we need to give the electromagnetic wave data on both sides. In addition, the thin layer problem can belong with multi-scale problems. The EM simulation is difficult to implement because of the high resource consumption of accurate modeling.</sec><sec>The accurate modeling is a very expensive algorithm, and the mesh-size difference between the thin layer region and the surrounding ordinary region is too large, thus seriously affecting the quality of the elements (tetrahedrons). So, it is difficult to solve the problem of thin layer quickly by using the time-domain discontinuous Galerkin method or even by using local time steps or other optimization methods. To acquire a fast solution, the approximate boundary condition can be used to describe the thin layer. The DGTD has also the advantage of the geometric fitting property. If a thin layer model is curved, we can build a conformal geometric model to deal with the thin layer, thus avoiding the complicated interpolation correction scheme.</sec><sec>In this work, A robust DGTD method for penetrable thin layers is proposed. In our study, we find that the current excitation scheme is prone to divergence when large dielectric coefficients are dealt with and its robustness is poor. To solve this problem, we use the electromagnetic field correction scheme in DGTD to deal with the thin layer problems. In addition, engineering problems may encounter multiple thin layers. In this work, a simulation scheme for multiple thin layers is presented and relevant numerical examples are given to verify it. The results illustrate the accuracy and effectiveness of the method in this work, and the method presented in this work can be used for the fast calculation of complex thin layer problems.</sec>

Novel Parallelization of Discontinuous Galerkin Method for Transient Electromagnetics Simulation Based on Sunway Supercomputers

A novel parallelization of discontinuous Galerkin time-domain (DGTD) method hybrid with the local time step (LTS) method on Sunway supercomputers for electromagnetic simulation is proposed. The proposed method includes a minimum number of roundtrip (MNR) strategy for processor-level parallelism and a double buffer strategy based on the remote memory access (RMA) of the Sunway processor. The MNR strategy optimizes the communication topology between nodes by recursively establishing the minimum spanning tree and the double buffer strategy is designed to make communication overlapped computation when RMA transmission. Combining the two methods, the proposed method achieves an unprecedented massively parallelism of the DGTD method. Several examples of radiation and scattering are used as cases to study the accuracy and validity of the proposed method. The numerical results show that the proposed method can effectively support 16,000 nodes (1,040,000 cores) parallelism on the Sunway supercomputer, which enables the DGTD method to solve the transient electromagnetic field in a very short time.

Uncertainty Quantification by Multilevel Monte Carlo and Local Time-Stepping for Wave Propagation

Uncertainty Quantification by Multilevel Monte Carlo and Local Time-Stepping for Wave Propagation

H-AMR: A New GPU-accelerated GRMHD Code for Exascale Computing with 3D Adaptive Mesh Refinement and Local Adaptive Time Stepping

General relativistic magnetohydrodynamic (GRMHD) simulations have revolutionized our understanding of black hole accretion. Here, we present a GPU-accelerated GRMHD code H-AMR with multifaceted optimizations that, collectively, accelerate computation by 2–5 orders of magnitude for a wide range of applications. First, it introduces a spherical grid with 3D adaptive mesh refinement that operates in each of the three dimensions independently. This allows us to circumvent the Courant condition near the polar singularity, which otherwise cripples high-resolution computational performance. Second, we demonstrate that local adaptive time stepping on a logarithmic spherical-polar grid accelerates computation by a factor of ≲10 compared to traditional hierarchical time-stepping approaches. Jointly, these unique features lead to an effective speed of ∼109 zone cycles per second per node on 5400 NVIDIA V100 GPUs (i.e., 900 nodes of the OLCF Summit supercomputer). We illustrate H-AMR's computational performance by presenting the first GRMHD simulation of a tilted thin accretion disk threaded by a toroidal magnetic field around a rapidly spinning black hole. With an effective resolution of 13,440 × 4608 × 8092 cells and a total of ≲22 billion cells and ∼0.65 × 108 time steps, it is among the largest astrophysical simulations ever performed. We find that frame dragging by the black hole tears up the disk into two independently precessing subdisks. The innermost subdisk rotation axis intermittently aligns with the black hole spin, demonstrating for the first time that such long-sought alignment is possible in the absence of large-scale poloidal magnetic fields.

On the Effect of Grid Orthogonalization in Stability and Accuracy of an FDTD Subgridding Method

We present a general-purpose subgridding (SG) inspired on a simplified orthogonalized integral (OI)-based method, leading to a simple and efficient algorithm with strong stability. It is based on a 2:1 transition ratio between coarse and finer zones and a natural local time-stepping (LTS) strategy to connect the domains. We provide a closed form for the stability condition and prove it heuristically, to show that roughly a 65% time-step reduction is enough to achieve stability. Results for coupling through lossy thin shells, compared with classical subcell impedance boundary methods, confirm the high accuracy of this method.

Performance Analysis of Speculative Parallel Adaptive Local Timestepping for Conservation Laws

Stable simulation of conservation laws, such as those used to model fluid dynamics and plasma physics applications, requires the satisfaction of the so-called Courant-Friedrichs-Lewy condition. By allowing regions of the mesh to advance with different timesteps that locally satisfy this stability constraint, significant work reduction can be attained when compared to a time integration scheme using a single timestep size. However, parallelizing this algorithm presents considerable difficulty. Since the stability condition depends on the state of the system, dependencies become dynamic and potentially non-local. In this article, we present an adaptive local timestepping algorithm using an optimistic (Timewarp-based) parallel discrete event simulation. We introduce waiting heuristics to limit misspeculation and a semi-static load balancing scheme to eliminate load imbalance as parts of the mesh require finer or coarser timesteps. Last, we outline an interface for separating the physics of the specific conservation law from the temporal integration allowing for productive adoption of our proposed algorithm. We present a misspeculation study for three conservation laws, demonstrating both the productivity of the local timestepping API, for which 74% of the lines of code are reused across different conservation laws, and the robustness of the waiting heuristics—at most 1.5% of element updates are rolled back. Our performance studies demonstrate up to a 2.8× speedup versus a baseline unoptimized local timestepping approach, a 4x improvement in per-node throughput compared to an MPI parallelization of synchronous timestepping, and scalability up to 3,072 cores on NERSC’s Cori Haswell partition.

An asynchronous variational integrator for the phase field approach to dynamic fracture

AbstractThe phase field approach is widely used to model fracture behaviors due to the absence of the need to track the crack topology and the ability to predict crack nucleation and branching. In this work, the asynchronous variational integrator (AVI) is adapted for the phase field approach of dynamic brittle fracture. The AVI is derived from Hamilton's principle and allows each element in the mesh to have its own local time step that may be different from others'. While the displacement field is explicitly updated, the phase field is implicitly solved, with upper and lower bounds strictly and conveniently enforced. In particular, two important variants of the phase field approach, the AT1 and AT2 models, are equally easily implemented. Several benchmark problems are used to study the performances of both the AT1 and AT2 models, and the results show that the AVI for the phase field approach significantly speeds up the computational efficiency and successfully captures the complicated dynamic fracture behavior.

Численное моделирование процессов дросселирования прямоточного воздушно-реактивного двигателя

Aimed at determining the maximum allowable pressure in the combustion chamber causing no flow separation in the air intake device, numerical simulation of the throttling process in a ramjet engine was carried out. Calculations were made in two formulations: with local and global time steps. When comparing results obtained using these approaches, the possibility was found to determine the air intake device throttling characteristics without losing the accuracy at the initial stages of designing the ramjet engines with introduction of the local time step. The non-stationary model created the basis for selecting the optimal integration step to minimize computational resources. The flow pulsations period was evaluated for the event of a surge. These pulsations frequencies were registered in the extremely low frequency range.

A novel transient-adaptive subcell algorithm with a hybrid application of different collision techniques in direct simulation Monte Carlo (DSMC)

A novel hybrid transient adaptive subcell (TAS) direct simulation Monte Carlo (DSMC) algorithm is proposed to simulate rarefied gas flows in a wide range of Knudsen numbers. It is derived and analyzed by using a time and spatial discrete operator approach based on the non-homogeneous, local N-particle kinetic equation, first proposed by Stefanov. The novel algorithm is considered together with the standard and hybrid collision algorithms built on uniform grids. The standard collision algorithm uses only one single scheme—the NoTime Counter (NTC), or the Generalized or Simplified Bernoulli trials (GBT, SBT). The hybrid algorithm employs NTC, GBT, or SBT depending on the instantaneous number of particles in the considered cell. The novel hybrid TAS algorithm benefits from both the hybrid collision approach and the transient adaptive subcell grid covering each collision cell to achieve a uniform accuracy of order O(Δt, Δr) independently of the number of particles in the cells. To this aim, a local time step is defined as coherent with the TAS grid covering the corresponding collision cell. The novel hybrid TAS algorithm is tested on two-dimensional benchmark problems: supersonic rarefied gas flow past of a flat plate under an angle of incidence and pressure-driven gas flow in a microchannel. The results obtained by the hybrid TAS algorithm are compared to those obtained by the standard algorithms and the available Bird's DS2V code using nearest neighbor collision and open-source OpenFOAM code. The comparison shows an excellent accuracy of the suggested algorithm in predicting the flow field.

Discontinuous Galerkin FEM with Hot Element Addition for the Thermal Simulation of Additive Manufacturing

Despite its promising advantages, the application of directed energy deposition (DED) to produce large metal parts is hindered by challenges inherent to the process. Undesired residual stresses, distortions and heterogeneous material properties mainly originate from a part’s thermal history. Fast part-scale thermal models therefore facilitate improved applicability of DED by enabling the prediction and mitigation of these unwanted effects. In this work, the efficiency of a discontinuous Galerkin-based thermal model with heat input by hot element addition, is evaluated and improved to allow such fast simulations. It is found that the model permits the use of a coarse discretization around the heat source, which significantly reduces simulation time while maintaining accurate solutions. It is also shown that the model naturally facilitates the use of local time stepping, which can considerably improve the efficiency of thermal additive manufacturing simulations.

A new two-phase shallow water hydro-sediment-morphodynamic model based on the HLLC solver and the hybrid LTS/GMaTS approach

In this paper, we present a new computationally-efficient and high-resolution depth-averaged two-phase flow model for hydro-sediment-morphydynamic processes, featuring an advance over existing models in terms of accuracy and efficiency of numerical solution. Under the framework of finite volume method (FVM) on unstructured grids, the Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver is proposed to compute inter-cell fluxes by applying the classical upwind HLLC approach to the water-sediment mixture and the sediment phases separately, in contrast to previous two-phase flow models using centered schemes. Moreover, to improve computational efficiency, the local-time-stepping (LTS) approach is implemented, the first attempt in the field of two-phase flow modelling. After a convergence rate study, the model is tested against a series of flow-sediment-bed evolutions induced by (1) two refilling processes of dredged trenches, (2) two instantaneous dam-break flooding flows, and (3) one levee breaching process by overtopping flows. It features encouraging performance when compared to a two-phase flow model based on a centered scheme and global time bound, characterized by more accurate results and much less computational cost. The present modelling framework shows promise in practical shallow water hydro-sediment-morphodynamic modelling applications.

Adaptive conservative time integration for transport in fractured porous media

The time step of an explicit time integration scheme for solving time-dependent hyperbolic partial differential equations in a finite volume method (FVM) framework is restricted by the Courant–Friedrichs–Lewy (CFL) criterion. Conventional time stepping integrates all grid cells with the same time step. This causes unnecessary computational costs when wave speeds and/or grid spacing vary considerably throughout the domain, and a few critical cells dictate the global time step, although most cells could be advanced with a much larger time step. Adaptive time stepping overcomes this issue by allowing different local time step sizes for each grid cell. The adaptive conservative time integration (ACTI) scheme is a recently developed adaptive time stepping method which relies on local time steps that are equal to the largest time step divided by powers of two.This work extends the ACTI scheme to tracer transport in fractured porous media. When fluid velocity within highly permeable fractures is higher than in the rock matrix and local grid refinement is applied around fractures, the CFL restriction would require prohibitively small time steps in the vicinity of fractures. For two-dimensional discrete fracture and matrix models of fracture patterns we demonstrate that ACTI reduces the computational cost by orders of magnitude compared to global time stepping while retaining solution accuracy. Empirically, we show that ACTI is stable in combination with a first-order explicit flux discretization scheme. Since combination with a standard higher-order MUSCL scheme can lead to spurious oscillations in the solution, we propose a modified MUSCL scheme relying on advection of an inclined reconstruction (MUSCL-AIR).

Advanced Parallelism of DGTD Method With Local Time Stepping Based on Novel MPI + MPI Unified Parallel Algorithm

In this communication, a novel message passing interface (MPI) parallel algorithm for nodal discontinuous Galerkin time-domain (NDGTD) method has been developed. A unified MPI + MPI technique has been introduced for extreme parallelism on a large-scale computer cluster. Through the data transmission between CPU nodes using MPI persistent nonblocking two-side communication and the direct data connection between processors in the same node via MPI shared memory windows, a two-layered parallel architecture is implemented to minimize the communication. To further accelerate the solution of the multiscale problems, the local time stepping (LTS) technique has been employed in the NDGTD method. A fast time step estimation method has been presented in this communication. With high overlap between the information transmission and the data calculation, the proposed MPI + MPI scheme overcomes the degradation of the parallel efficiency of the pure MPI technique in the scenario of the LTS technique and the large-scale CPU cores. Up to 94% parallel efficiency in 6400 CPU cores is achieved for the average single-core loading about 1700 finite elements, and 18 times acceleration for time step estimation can be obtained with the fourth-order basis function. Three practical complex examples are given to demonstrate a good performance of the proposed method.

Adaptive total variation stable local timestepping for conservation laws

This paper proposes a first-order total variation diminishing (TVD) treatment for coarsening and refining of local timestep size in response to dynamic local variations in wave speeds for nonlinear conservation laws. The algorithm is accompanied with a proof of formal correctness showing that given a sufficiently small minimum timestep the algorithm will produce a TVD solution for nonlinear scalar conservation laws. A key feature of the algorithm is its formulation as a discrete event simulation, which allows for easy and efficient parallelization using existing software. Numerical results demonstrate the stability and adaptivity of the method for the shallow water equations. We also introduce a performance model to load balance and explain the observed performance gains. Performance results are presented for a single node on Stampede2's Skylake partition using an optimistic parallel discrete event simulator. Results show the proposed algorithm recovering 59%-77% of the theoretically achievable speed-up with the remainder being attributed to the cost of computing the CFL condition and load imbalance.

Scalable Local Timestepping on Octree Grids

Numerical solutions of hyperbolic partial differential equations (PDEs) are ubiquitous in science and engineering. Method of lines is a popular approach to discretize PDEs defined in spacetime, where space and time are discretized independently. When using explicit timesteppers on adaptive grids, the use of a global timestep size dictated by the finest grid-spacing leads to inefficiencies in the coarser regions. Even though adaptive space discretizations are widely used in computational sciences, temporal adaptivity is less common due to its sophisticated nature. This paper presents highly scalable algorithms to enable local timestepping (LTS) for explicit timestepping schemes on octrees. We demonstrate the accuracy of our methods and the scalability of our framework across 16K cores on TACC's Frontera supercomputer https://frontera-portal.tacc.utexas.edu/. We also present a speedup estimation model for LTS that predicts the distributed speedup compared to global timestepping with an average of 0.1 relative error.

A GPU-Accelerated and LTS-Based Finite Volume Shallow Water Model

This paper presents a GPU (Graphics Processing Unit)-accelerated and LTS (Local-time-Step)-based finite volume Shallow Water Model (SWM). The model performance is compared against the other five model versions (Single CPU versions with/without LTS, Multi-CPU versions with/without LTS, and a GPU version) by simulating three flow scenarios: an idealized dam-break flow; an experimental dam-break flow; a field-scale scenario of tidal flows. Satisfactory agreements between simulation results and the available measured data/reference solutions (water level, flow velocity) indicate that all the six SWM versions can well simulate these challenging shallow water flows. Inter-comparisons of the computational efficiency of the six SWM versions indicate the following. First, GPU acceleration is much more efficient than multi-core CPU parallel computing. Specifically, the speed increase in the GPU can be as high as a hundred, whereas those for multi-core CPU are only 2–3. Second, implementing the LTS can bring considerable reduction: the additional maximum speed-ups can be as high as 10 for the single-core CPU/multi-core CPU versions, and as high as five for the GPU versions. Third, the GPU + LTS version is computationally the most efficient in most cases; the multi-core CPU + LTS version may run as fast as a GPU version for scenarios over some intermediate number of cells.

A parallel modular computing environment for three-dimensional multiresolution simulations of compressible flows

Numerical investigation of compressible flows faces two main challenges. In order to accurately describe the flow characteristics, high-resolution nonlinear numerical schemes are needed to capture discontinuities and resolve wide convective, acoustic and interfacial scale ranges. The simulation of realistic three-dimensional (3D) problems with state-of-the-art finite-volume method (FVM) based on approximate Riemann solvers with weighted nonlinear reconstruction schemes requires the usage of high-performance computing (HPC) architectures. Efficient compression algorithms reduce computational and memory load. Fully adaptive multiresolution (MR) algorithms with LTS have proven their potential for such applications. While modern central processing units (CPUs) requires multiple levels of parallelism to achieve peak performance, the fine-grained MR mesh adaptivity results in challenging compute/communication patterns. Moreover, local time stepping (LTS) incurs strong data dependencies which challenge a parallelization strategy.We address these challenges with a block-based MR algorithm, where arbitrary cuts in the underlying octree are possible. This allows for a parallelization on distributed-memory machines via the Message Passing Interface (MPI). We obtain neighbor relations by a simple bit-logic in a modified Morton Order. The block-based concept allows for a modular setup of the source code framework in which the building blocks of the algorithm, such as the choice of the Riemann solver or the reconstruction stencil, are interchangeable without loss of parallel performance. We present the capabilities of the modular framework with a range of test cases and scaling analysis with effective resolutions beyond one billion cells using O(104) cores.