Implicit Time Stepping Research Articles

On the convergence of a linearly implicit finite element method for the nonlinear Schrödinger equation

AbstractWe consider a model initial‐ and Dirichlet boundary–value problem for a nonlinear Schrödinger equation in two and three space dimensions. The solution to the problem is approximated by a conservative numerical method consisting of a standard conforming finite element space discretization and a second‐order, linearly implicit time stepping, yielding approximations at the nodes and at the midpoints of a nonuniform partition of the time interval. We investigate the convergence of the method by deriving optimal‐order error estimates in the and the norm, under certain assumptions on the partition of the time interval and avoiding the enforcement of a Courant‐Friedrichs‐Lewy (CFL) condition between the space mesh size and the time step sizes.

Approximation‐based implicit integration algorithm for the Simo‐Miehe model of finite‐strain inelasticity

AbstractWe propose a simple, efficient, and reliable procedure for implicit time stepping, regarding a special case of the viscoplasticity model proposed by Simo and Miehe (1992). The kinematics of this popular model is based on the multiplicative decomposition of the deformation gradient tensor, allowing for a combination of Newtonian viscosity and arbitrary isotropic hyperelasticity. The algorithm is based on approximation of precomputed solutions. Both Lagrangian and Eulerian versions of the algorithm with equivalent properties are available. The proposed numerical scheme is non‐iterative, unconditionally stable, and first order accurate. Moreover, the integration algorithm strictly preserves the inelastic incompressibility constraint, symmetry, positive definiteness, and w‐invariance. The accuracy of stress calculations is verified in a series of numerical tests, including non‐proportional loading and large strain increments. In terms of stress calculation accuracy, the proposed algorithm is equivalent to the implicit Euler method with strict inelastic incompressibility. The algorithm is implemented into MSC.MARC and a demonstration initial‐boundary value problem is solved.

Accelerating aeroelastic UVLM simulations by inexact Newton algorithms

AbstractWe consider the aeroelastic simulation of flexible mechanical structures submerged in subsonic fluid flows at low Mach numbers. The nonlinear kinematics of flexible bodies are described in the total Lagrangian formulation and discretized by finite elements. The aerodynamic loads are computed using the unsteady vortex-lattice method wherein a free wake is tracked over time. Each implicit time step in the dynamic simulation then requires solving a nonlinear equation system in the structural variables with additional aerodynamic load terms. Our focus here is on the efficient numerical solution of this system by accelerating the Newton algorithm. The particular structure of the aeroelastic nonlinear system suggests the structural derivative as an approximation to the full derivative in the linear Newton system. We investigate and compare two promising algorithms based on this approximation, a quasi-Newton type algorithm and a novel inexact Newton algorithm. Numerical experiments are performed on a flexible plate and on a wind turbine. Our computational results show that the approximation can indeed accelerate the Newton algorithm substantially. Surprisingly, the theoretically preferable inexact Newton algorithm is much slower than the quasi-Newton algorithm, which motivates further research to speed up derivative evaluations.

ReMKiT1D - A framework for building reactive multi-fluid models of the tokamak scrape-off layer with coupled electron kinetics in 1D

In this manuscript we present the recently developed flexible framework for building both fluid and electron kinetic models of the tokamak Scrape-Off Layer in 1D - ReMKiT1D (Reactive Multi-fluid and Kinetic Transport in 1D). The framework can handle systems of non-linear ODEs, various 1D PDEs arising in fluid modelling, as well as PDEs arising from the treatment of the electron kinetic equation. As such, the framework allows for flexibility in fluid models of the Scrape-Off Layer while allowing the easy addition of kinetic electron effects. We focus on presenting both the high-level design decisions that allow for model flexibility, as well as the most important implementation aspects. A significant number of verification and performance tests are presented, as well as a step-by-step walkthrough of a simple example for setting up models using the Python interface. Program summaryProgram title: ReMKiT1DCPC Library link to program files:https://doi.org/10.17632/j47ym66xzj.1Developer's repository link:https://github.com/ukaea/ReMKiT1D and https://github.com/ukaea/ReMKiT1D-PythonLicensing provisions: GPLv3Programming language: Fortran, PythonSupplementary material:https://doi.org/10.14468/fdq7-z869Nature of problem: The flexible generation and modification of 1D models pertaining to multi-fluid simulations of the tokamak Scrape-Off Layer (SOL) with electron kinetics and reaction support. This would then allow both for rapid iteration on reduced models as well as the evaluation of kinetic electron effects in equilibria and during transients, following the formalism previously developed for SOL-KiT [1]. The framework was not only envisioned as the successor to SOL-KiT, but a tool that would allow users to construct their own models coupled with electron kinetics capabilities.Solution method: The framework is written heavily utilizing Object-Oriented design principles, in particular using an extended version of the puppeteer pattern as presented by Rouson et al. [2], as well as the heavy use of the strategy pattern/dependency injection. The Fortran code is MPI parallel and utilizes the PETSc library for implicit time-stepping. MPI parallelization is extended to distribution function Legendre harmonics, allowing for improved strong scaling. Initialization of the Fortran framework is done using JSON configuration files generated by an accompanying Python interface, and data analysis is standardized using widely used data formats such as HDF5.Additional comments including restrictions and unusual features: The present manuscript focuses on the design and high-level implementation of the framework, as well as the demonstration of the workflow and various verification and performance benchmarking tests. Some details are avoided for the sake of brevity at various points, and these are meant to be available as part of the general code documentation or tutorials offered on the main repositories.

Scalable implicit solvers with dynamic mesh adaptation for a relativistic drift-kinetic Fokker–Planck–Boltzmann model

In this work we consider a relativistic drift-kinetic model for runaway electrons along with a Fokker–Planck operator for small-angle Coulomb collisions, a radiation damping operator, and a secondary knock-on (Boltzmann) collision source. We develop a new scalable fully implicit solver utilizing finite volume and conservative finite difference schemes and dynamic mesh adaptivity. A new data management framework in the PETSc library based on the p4est library is developed to enable simulations with dynamic adaptive mesh refinement (AMR), distributed memory parallelization, and dynamic load balancing of computational work. This framework and the runaway electron solver building on the framework are able to dynamically capture both bulk Maxwellian at the low-energy region and a runaway tail at the high-energy region. To effectively capture features via the AMR algorithm, a new AMR indicator prediction strategy is proposed that is performed alongside the implicit time evolution of the solution. This strategy is complemented by the introduction of computationally cheap feature-based AMR indicators that are analyzed theoretically. Numerical results quantify the advantages of the prediction strategy in better capturing features compared with nonpredictive strategies; and we demonstrate trade-offs regarding computational costs. The robustness with respect to model parameters, algorithmic scalability, and parallel scalability are demonstrated through several benchmark problems including manufactured solutions and solutions of different physics models. We focus on demonstrating the advantages of using implicit time stepping and AMR for runaway electron simulations.

Accuracy Analysis for Explicit-Implicit Finite Volume Schemes on Cut Cell Meshes

The solution of time-dependent hyperbolic conservation laws on cut cell meshes causes the small cell problem: standard schemes are not stable on the arbitrarily small cut cells if an explicit time stepping scheme is used and the time step size is chosen based on the size of the background cells. In May and Berger (J Sci Comput 71: 919–943, 2017), the mixed explicit-implicit approach in general and MUSCL-Trap (monotonic upwind scheme for conservation laws and trapezoidal scheme) in particular have been introduced to solve this problem by using implicit time stepping on the cut cells. Theoretical and numerical results have indicated that this might lead to a loss in accuracy when switching between the explicit and implicit time stepping. In this contribution, we examine this in more detail and will prove in one dimension that the specific combination MUSCL-Trap of an explicit second-order and an implicit second-order scheme results in a fully second-order mixed scheme. As this result is unlikely to hold in two dimensions, we also introduce two new versions of mixed explicit-implicit schemes based on exchanging the explicit scheme. We present numerical tests in two dimensions where we compare the new versions with the original MUSCL-Trap scheme.

Two arbitrary-order constraint-preserving schemes for the Yang–Mills equations on polyhedral meshes

<abstract><p>Two numerical schemes are proposed and investigated for the Yang–Mills equations, which can be seen as a nonlinear generalisation of the Maxwell equations set on Lie algebra-valued functions, with similarities to certain formulations of General Relativity. Both schemes are built on the Discrete de Rham (DDR) method, and inherit from its main features: an arbitrary order of accuracy, and applicability to generic polyhedral meshes. They make use of the complex property of the DDR, together with a Lagrange-multiplier approach, to preserve, at the discrete level, a nonlinear constraint associated with the Yang–Mills equations. We also show that the schemes satisfy a discrete energy dissipation (the dissipation coming solely from the implicit time stepping). Issues around the practical implementations of the schemes are discussed; in particular, the assembly of the local contributions in a way that minimises the price we pay in dealing with nonlinear terms, in conjunction with the tensorisation coming from the Lie algebra. Numerical tests are provided using a manufactured solution, and show that both schemes display a convergence in $ L^2 $-norm of the potential and electrical fields in $ \mathcal O(h^{k+1}) $ (provided that the time step is of that order), where $ k $ is the polynomial degree chosen for the DDR complex. We also numerically demonstrate the preservation of the constraint.</p></abstract>

Robust conservative scheme and nonlinear solver for phase transitions in heterogeneous permafrost

In this paper we study computational schemes to simulate freezing and thawing in permafrost modeled by a nonlinear heat equation with constitutive properties resembling those in the well known Stefan problem but featuring distinct challenges. The models are discretized with low order conservative discretization and fully implicit time stepping. We explore the challenges due to the nonsmooth nature of the temperature–enthalpy relationship and the domain heterogeneity, with focus on the solver and the use of enthalpy as a primary variable in contrast to the temperature variable used commonly in applications literature. We prove and demonstrate the convergence of our algorithms in realistic physical scenarios.

A compatible finite element discretisation for the nonhydrostatic vertical slice equations

We present a compatible finite element discretisation for the vertical slice compressible Euler equations, at next-to-lowest order (i.e., the pressure space is bilinear discontinuous functions). The equations are numerically integrated in time using a fully implicit timestepping scheme which is solved using monolithic GMRES preconditioned by a linesmoother. The linesmoother only involves local operations and is thus suitable for domain decomposition in parallel. It allows for arbitrarily large timesteps but with iteration counts scaling linearly with Courant number in the limit of large Courant number. This solver approach is implemented using Firedrake, and the additive Schwarz preconditioner framework of PETSc. We demonstrate the robustness of the scheme using a standard set of testcases that may be compared with other approaches.

A Linear and Angular Momentum Conserving Hybrid Particle/Grid Iteration for Volumetric Elastic Contact

We present a momentum conserving hybrid particle/grid iteration for resolving volumetric elastic collision. Our hybrid method uses implicit time stepping with a Lagrangian finite element discretization of the volumetric elastic material together with impulse-based collision-correcting momentum updates designed to exactly conserve linear and angular momentum. We use a two-step process for collisions: first we use a novel grid-based approach that leverages the favorable collision resolution properties of Particle-In-Cell (PIC) techniques, then we finalize with a classical collision impulse strategy utilizing continuous collision detection. Our PIC approach uses Affine-Particle-In-Cell momentum transfers as collision preventing impulses together with novel perfectly momentum conserving boundary resampling and downsampling operators that prevent artifacts in portions of the boundary where the grid resolution is of disparate resolution. We combine this with a momentum conserving augury iteration to remove numerical cohesion and model sliding friction. Our collision strategy has the same continuous collision detection as traditional approaches, however our hybrid particle/grid iteration drastically reduces the number of iterations required. Lastly, we develop a novel symmetric positive semi-definite Rayleigh damping model that increases the convexity of the nonlinear systems associated with implicit time stepping. We demonstrate the robustness and efficiency of our approach in a number of collision intensive examples.

High-order variational Lagrangian schemes for compressible fluids

We present high-order variational Lagrangian finite element methods for compressible fluids using a discrete energetic variational approach. Our spatial discretization is mass/momentum/energy conserving and entropy stable. Fully implicit time stepping is used for the temporal discretization, which allows for a much larger time step size for stability compared to explicit methods, especially for low-Mach number flows and/or on highly distorted meshes. Ample numerical results are presented to showcase the good performance of our proposed scheme.

Time parallelism and Newton-adaptivity of the two-derivative deferred correction discontinuous Galerkin method

In this work, we consider a high-order discretization of compressible viscous flows allowing parallelization both in space and time.The discontinuous Galerkin spectral element method, which is well-suited for massively parallel simulations, is used for spatial discretization. The main novelty in this work is the additional demonstration of time-parallel capabilities within an implicit two-derivative timestepping procedure to further increase the parallel speedup. Temporal parallelism is made possible by a predictor-corrector-type time discretization that allows to split the associated workload onto multiple processors.We identify a homogeneous load balance with respect to the linear (GMRES) iterations on each processor as a key for parallel efficiency. To homogenize the load and to enable practical simulations, an adaptive strategy for Newton’s method is introduced. It is shown that the time-parallel method provides a parallel efficiency of approx. 60−70% on 4−7 computational partitions. Moreover, the capabilities of the novel method for the simulation of large-scale problems are illustrated with a mixed temporal and spatial parallelization on more than 1000 processors.

Extensional flow of a free film of nematic liquid crystal with moderate elasticity

The human tear film is a multilayer structure in which the dynamics are often strongly affected by a floating lipid layer. That layer has liquid crystalline characteristics and plays important roles in the health of the tear film. Previous models have treated the lipid layer as a Newtonian fluid in extensional flow. Motivated to develop a more realistic treatment, we present a model for the extensional flow of thin sheets of nematic liquid crystal. The rod-like molecules of these substances impart an elastic contribution to the rheology. We rescale a weakly elastic model due to Cummings et al. [“Extensional flow of nematic liquid crystal with an applied electric field,” Eur. J. Appl. Math. 25, 397–423 (2014).] to describe a lipid layer of moderate elasticity. The resulting system of two nonlinear partial differential equations for sheet thickness and axial velocity is fourth order in space, but still represents a significant reduction of the full system. We analyze solutions arising from several different boundary conditions, motivated by the underlying application, with particular focus on dynamics and underlying mechanisms under stretching. We solve the system numerically, via collocation with either finite difference or Chebyshev spectral discretization in space, together with implicit time stepping. At early times, depending on the initial film shape, pressure either aids or opposes extensional flow, which changes the free surface dynamics of the sheet and can lead to patterns reminiscent of those observed in tear films. We contrast this finding with the cases of weak elasticity and Newtonian flow, where the sheet retains the same qualitative shape throughout time.

Adaptively implicit MPDATA advection for arbitrary Courant numbers and meshes

AbstractAdvection schemes with time step restrictions are widely used in weather and climate models. This can lead to instability in the presence of high flow speeds (relative to mesh spacing) such as occurs in convective updraughts, regions of mesh convergence, or where the winds are unusually high. An adaptively implicit advection scheme is proposed that treats advection implicitly only where the Courant number is high. Flux correction to ensure monotonicity is adapted to work with implicit time stepping. A version of the multidimensional positive‐definite advection transport algorithm MPDATA is derived with an anti‐diffusive flux compensating truncation errors of off‐centred implicit time stepping. The anti‐diffusive flux is gradually reduced as Courant numbers increase above 2 in order to maintain stability at the expense of second‐order accuracy at high Courant numbers. Results of two‐dimensional advection by deformational flow are presented on various meshes of the sphere. Stability and second‐order accuracy are maintained when the Courant number is over 100 in a small region, when strong wind crosses the poles of a rotated latitude–longitude mesh. Good solutions are also obtained on a skipped latitude–longitude mesh, a cubed sphere, and hexagonal meshes. Accuracy reverts to first order when Courant numbers are large over a large fraction of the domain.

Discrete gradient flow approximations of high dimensional evolution partial differential equations via deep neural networks

We consider the approximation of initial/boundary value problems involving, possibly high-dimensional, dissipative evolution partial differential equations (PDEs) using a deep neural network framework. More specifically, we first propose discrete gradient flow approximations based on non-standard Dirichlet energies for problems involving essential boundary conditions posed on bounded spatial domains. The imposition of the boundary conditions is realized weakly via non-standard functionals; the latter classically arise in the construction of Galerkin-type numerical methods and are often referred to as “Nitsche-type” methods. Moreover, inspired by the seminal work of Jordan, Kinderlehrer, and Otto (JKO) Jordan et al. (1998), we consider the second class of discrete gradient flows for special classes of dissipative evolution PDE problems with non-essential boundary conditions. These JKO-type gradient flows are solved via deep neural network approximations. A key, distinct aspect of the proposed methods is that the discretization is constructed via a sequence of residual-type deep neural networks (DNN) corresponding to implicit time-stepping. As a result, a DNN represents the PDE problem solution at each time node. This approach offers several advantages in the training of each DNN. We present a series of numerical experiments which showcase the good performance of Dirichlet-type energy approximations for lower space dimensions and the excellent performance of the JKO-type energies for higher spatial dimensions.

Implicit ensemble tangent linear models (IETLMs) for model differentiation

AbstractIdeally, tangent linear models (TLMs) predict the difference between perturbed and unperturbed non‐linear forecasts of interest. The adjoint of a TLM gives the gradient of the non‐linear model and is used in 4DVar data assimilation and in adjoint‐based Forecast Sensitivity to Observation Impact (FSOI). The accuracy of the local ensemble TLM (LETLM) has been shown to be limited by its inability to account for implicit time stepping. Here we derive implicit ensemble TLMs (IETLMs) that, at most, require the number of independent ensemble members to be equal to the number of variables in the implicit computational stencil. The accuracy of the IETLM in the linear regime is confirmed using an implicitly time‐stepped Lorenz 96 model and a nine‐member ensemble. IETLMs feature two sparse matrices: matrix N that operates on an initially unknown future time perturbation, and matrix L that operates on the current time perturbation. For ensemble perturbations in the non‐linear regime, we develop a diagonally robust (DR) IETLM that reduces the chances of N becoming ill‐conditioned. The performance of the DR IETLM was compared with traditional TLM performance using IETLM ensemble perturbations whose ‘Gilmour et al., 2001’ measure of non‐linearity ranged up to the non‐linearity of operational 32‐hr ensemble forecast perturbations. Over a wide range of non‐linearity, the DR IETLM performance was found to match that of the traditional TLM provided the initial standard deviation of the ensemble perturbations was times the standard deviation of the test perturbations. Ideal FSOI requires the adjoint of a TLM that accurately predicts the known difference between corrected and uncorrected non‐linear forecasts. The DR IETLM was found to meet this FSOI accuracy requirement much more closely than the traditional TLM when the ensemble perturbations were created by subtracting the corrected forecast from ensemble members that were centred on the uncorrected forecast.

A two‐fluid single‐column model of turbulent shallow convection. Part II: Single‐column model formulation and numerics

AbstractThe two‐fluid single‐column model of Thuburn et al. (Quart. J. R. Meteorol. Soc., 2019, 145, 1535–1550) is extended to include moisture and horizontal wind shear. Turbulent kinetic energy is introduced as a prognostic variable, dependence on a diagnosed boundary‐layer height is removed, and subfilter fluxes are approximated using a two‐fluid version of a Mellor–Yamada scheme. Three mechanisms for entrainment and detrainment processes are introduced, which represent entrainment of unstable air at the surface, forced detrainment of air at the top of the boundary/cloud layers, and turbulent mixing that relaxes the convective fluid to a reference profile. A semi‐implicit Eulerian discretization replaces the semi‐implicit semi‐Lagrangian implementation of Thuburn et al. (Quart. J. R. Meteorol. Soc., 2019, 145, 1535–1550) to improve numerical stability and conservation. The equations for the implicit time step are solved using a quasi‐Newton method, which is shown to perform well in numerical tests for conservation and convergence. The two‐fluid single‐column model presented in this article will be applied to simulations of shallow cumulus convection in Part III.

Analysis of nonlinear soil-structure interaction using partitioned method

For the current state-of-practice in soil-structure interaction (SSI) of nuclear facilities, the strain-dependent characteristics of soil are considered indirectly via equivalent linear methods in the frequency domain, typically represented by the SASSI program. SSI analysis in the time-domain, although directly involving material nonlinearity of the soil, is inefficient in practice up to now. Based on Partitioned Analysis of SSI (PASSI), a method for nonlinear soil-structure interaction analysis is developed in this paper, by applying implicit time-step integration and explicit one with different time steps to nuclear island and to soil, respectively. The incremental equilibrium equation and explicit decoupling method are used to analyze the soil nonlinearity described by the Davidenkov model with simplified loading-reloading rules, which avoids solving algebraic equations and iteration processes. An asynchronous parallel algorithm for nonlinear SSI analysis is given. A simple example is given to verify the partitioned approach with explicit-implicit co-computation against the full-explicit approach. Seismic response characteristics of nuclear power plants are investigated by comparing responses of nonlinear analysis with linear analysis. Different results address the importance of nonlinear SSI analysis for the seismic performance of nuclear structures.

Parallel implementation of the SHYFEM (System of HydrodYnamic Finite Element Modules) model

Abstract. This paper presents the message passing interface (MPI)-based parallelization of the three-dimensional hydrodynamic model SHYFEM (System of HydrodYnamic Finite Element Modules). The original sequential version of the code was parallelized in order to reduce the execution time of high-resolution configurations using state-of-the-art high-performance computing (HPC) systems. A distributed memory approach was used, based on the MPI. Optimized numerical libraries were used to partition the unstructured grid (with a focus on load balancing) and to solve the sparse linear system of equations in parallel in the case of semi-to-fully implicit time stepping. The parallel implementation of the model was validated by comparing the outputs with those obtained from the sequential version. The performance assessment demonstrates a good level of scalability with a realistic configuration used as benchmark.

A Multiresolution Discrete Element Method for Triangulated Objects with Implicit Time stepping

Simulations of many rigid bodies colliding with each other sometimes yield particularly interesting results if the colliding objects differ significantly in size and are non-spherical. The most expensive part within such a simulation code is the collision detection. We propose a family of novel multiscale collision detection algorithms that can be applied to triangulated objects within explicit and implicit time stepping methods. They are well-suited to handle objects that cannot be represented by analytical shapes or assemblies of analytical objects. Inspired by multigrid methods and adaptive mesh refinement, we determine collision points iteratively over a resolution hierarchy, and combine a functional minimisation plus penalty parameters with the actual comparision-based geometric distance calculation. Coarse surrogate geometry representations identify "no collision" scenarios early on and otherwise yield an educated guess which triangle subsets of the next finer level potentially yield collisions. They prune the search tree, and furthermore feed conservative contact force estimates into the iterative solve behind an implicit time stepping. Implicit time stepping and non-analytical shapes often yield prohibitive high compute cost for rigid body simulations. Our approach reduces these cost algorithmically by one to two orders of magnitude. It also exhibits high vectorisation efficiency due to its iterative nature.