High-order Integration Research Articles

A robust upwind mixed hybrid finite element method for transport in variably saturated porous media

Abstract. The mixed finite element (MFE) method is well adapted for the simulation of fluid flow in heterogeneous porous media. However, when employed for the transport equation, it can generate solutions with strong unphysical oscillations because of the hyperbolic nature of advection. In this work, a robust upwind MFE scheme is proposed to avoid such unphysical oscillations. The new scheme is a combination of the upwind edge/face centered finite volume method with the hybrid formulation of the MFE method. The scheme ensures continuity of both advective and dispersive fluxes between adjacent elements and allows to maintain the time derivative continuous, which permits employment of high-order time integration methods via the method of lines (MOL). Numerical simulations are performed in both saturated and unsaturated porous media to investigate the robustness of the new upwind MFE scheme. Results show that, contrarily to the standard scheme, the upwind MFE method generates stable solutions without under and overshoots. The simulation of contaminant transport into a variably saturated porous medium highlights the robustness of the proposed upwind scheme when combined with the MOL for solving nonlinear problems.

Efficient adaptive step size control for exponential integrators

Traditional step size controllers make the tacit assumption that the cost of a time step is independent of the step size. This is reasonable for explicit and implicit integrators that use direct solvers. In the context of exponential integrators, however, an iterative approach, such as Krylov methods or polynomial interpolation, to compute the action of the required matrix functions is usually employed. In this case, the assumption of constant cost is not valid. This is, in particular, a problem for higher-order exponential integrators, which are able to take relatively large time steps based on accuracy considerations. In this paper, we consider an adaptive step size controller for exponential Rosenbrock methods that determines the step size based on the premise of minimizing computational cost. The largest allowed step size, given by accuracy considerations, merely acts as a constraint. We test this approach on a range of nonlinear partial differential equations. Our results show significant improvements (up to a factor of 4 reduction in the computational cost) over the traditional step size controller for a wide range of tolerances.

High-order integrators for Lagrangian systems on homogeneous spaces via nonholonomic mechanics

In this paper, high-order numerical integrators on homogeneous spaces will be presented as an application of nonholonomic partitioned Runge–Kutta Munthe-Kaas (RKMK) methods on Lie groups. A homogeneous space M is a manifold where a group G acts transitively. Such a space can be understood as a quotient M≅G/H, where a closed Lie subgroup H, is the isotropy group of each point of M. The Lie algebra of G may be decomposed into g=m⊕h, where h is the subalgebra that generates H and m is a subspace spanning the tangent spaces of M. Thus, variational problems on M can be treated as nonholonomically constrained problems on G, by requiring variations to remain on m. Nonholonomic partitioned RKMK integrators are derived as a modification of those obtained by a discrete variational principle on Lie groups, and can be interpreted as obeying a discrete Chetaev principle. These integrators tend to preserve several properties of their purely variational counterparts.

Arbitrary high-order exponential integrators conservative schemes for the nonlinear Gross-Pitaevskii equation

In this paper, we propose a family of high-order conservative schemes based on the exponential integrators technique and the symplectic Runge-Kutta method for solving the nonlinear Gross-Pitaevskii equation. By introducing generalized scalar auxiliary variable, the equation is equivalent to a new system with both mass and modified energy conservation laws. Then, a conservative semi-discrete exponential scheme is given by combining the symplectic Runge-Kutta method and the Lawson method in the time direction. Subsequently, we apply the Fourier pseudo-spectral method to approximate semi-discrete system in space and obtain the fully-discrete schemes that conserve the energy and mass. Numerical examples are presented to confirm the accuracy and conservation of the developed schemes.

Lagrangian coherent structures in flow past a backward-facing step

This paper investigates flow past a backward-facing step (BFS) in a duct at Reynolds number $Re = 5080$ based on step height, mean inflow velocity and fluid kinematic viscosity. The flow configuration matches a combustion experiment conducted by Pitz and Daily in 1983. High-resolution velocity fields are obtained in OpenFOAM by direct numerical simulation, and the flow field analysed by Lagrangian approaches. Trajectories of fluid particles in areas of interest are obtained by high-order numerical integration, and used to compute finite-time Lyapunov exponents (FTLEs) and polar rotation angles. Lagrangian coherent structures (LCSs) are extracted using geodesic theory, including hyperbolic LCSs and elliptic LCSs. We use complementary qualitative and quantitative LCS analyses to uncover the underlying flow structures. Notably, we find that a flow pathway in which fluid particles rarely diverge from adjacent particles is opened and closed by FTLE ridges determined by the periodic shedding of vortices from the BFS. Two dominant vortices with significant Lagrangian coherence, generated respectively by the separated boundary layer and shear layer, are self-sustaining and of comparable strength. Hyperbolic repelling LCSs act as transport barriers between the pathway and cores of the coherent vortices, thus playing a major part in the fluid entrainment process. Interactions between these different geometric regions partitioned by LCSs lead to intrinsic complexity in the BFS flow.

Exponential Bhatnagar-Gross-Krook integrator for multiscale particle-based kinetic simulations.

Despite the development of an extensive toolbox of multiscale rarefied flow simulators, such simulations remain challenging due to the significant disparity of collisional and macroscopic spatiotemporal scales. Our study offers a novel and consistent numerical scheme for a coupled treatment of particles advection and collision governed by the BGK evolution, honoring positivity of the velocity distribution. Our method shares its framework, in spirit, with the unified gas kinetic class of multiscale schemes. Yet it provides attractive features for particle-based stochastic simulations, readily implementable to existing direct simulation Monte Carlo codes. Two main innovations are integrated in the presented BGK particle method. The first ingredient is a high-order time integration that can be interpreted probabilistically, independent of the time step size. The next one is identifying modified particle distributions that remain invariant under the advection-relaxation evolution. We demonstrate accuracy and performance of the devised scheme for prototypic gas flows over a wide range of rarefaction parameters. Due to the resulting robustness and flexibility of the devised exponential BGK integrator, the scheme paves the way towards more affordable simulations of large-scale and multiscale rarefied gas phenomena.

An implicit SPH-based structure model for accurate Fluid–Structure Interaction simulations with hourglass control scheme

This work is dedicated to development of an implicit SPH (Smoothed Particle Hydrodynamics)-based structure model as well as its integrated purely Lagrangian meshfree hydroelastic FSI (Fluid–Structure Interaction) solver for consistent, accurate and stable modeling of FSI problems. The implicit structure model is established in the context of Hamiltonian mechanics under the assumption of linear elastic solid. Four refined schemes are introduced for improved accuracy, consistency, efficiency and stability: two formerly developed schemes, [i] HT (High-order implicit Time integration) scheme, [ii] HD (High-order Discretization) operator scheme; as well as the novel two schemes newly proposed in this work, namely [iii] IPC (Iterative Predictor–Corrector) calculation procedures to minimize errors related to an imprecise assumption on rotation matrix and [iv] IHC (Implicit Hourglass Control) scheme for stabilizing simulation by suppressing spurious zero-energy modes. An enhanced particle-based FSI solver is formulated within a unified SPH framework by integrating a refined ISPH (Incompressible SPH) fluid model and the proposed implicit structure model, which allows consistent fluid–structure time coupling adopting the equivalent time step sizes in both phases. The proposed implicit structure model and FSI solver are configured in two-dimensions and validated with several 2D structure and FSI benchmark tests showing that our proposed implicit framework could provide almost consistent robustness, accuracy and efficiency with respect to the explicit one.

Conservative integrators for vortex blob methods on the plane

Conservative integrators have recently been derived for many-body systems using the Discrete Multiplier Method, specifically for the point vortex problems in [1]. In the present paper, we extend this work to derive conservative symmetric second-order one-step integrators for a family of vortex blob models. A rational function approximation was used to approximate the exponential integral which appears in the Hamiltonian. Conservative properties and second order convergence are proved. Numerical experiments are shown to verify the conservative property of these integrators, their second-order accuracy, and as well as the resulting spatial and temporal accuracy of the vortex blob method. Moreover, the derived implicit conservative integrators are shown to better preserve conserved quantities than standard higher-order explicit integrators on comparable computation times.

Personalized Algorithm Generation: A Case Study in Learning ODE Integrators

We study the learning of numerical algorithms for scientific computing, which combines mathematically driven, handcrafted design of general algorithm structure with a data-driven adaptation to specific classes of tasks. This represents a departure from the classical approaches in numerical analysis, which typically do not feature such learning-based adaptations. As a case study, we develop a machine learning approach that automatically learns effective solvers for initial value problems in the form of ordinary differential equations (ODEs), based on the Runge--Kutta (RK) integrator architecture. We show that we can learn high-order integrators for targeted families of differential equations without the need for computing integrator coefficients by hand. Moreover, we demonstrate that in certain cases we can obtain superior performance to classical RK methods. This can be attributed to certain properties of the ODE families being identified and exploited by the approach. Overall, this work demonstrates an effective learning-based approach to the design of algorithms for the numerical solution of differential equations. This can be readily extended to other numerical tasks.

FREQUENCY FEATURES OF THE NUMERICAL METHOD SAMPLING OF DIGITAL CONTROL SYSTEMS

Context. The studies of the frequency properties of the explicit multistep numerical integrators which use for sampling of continuous transfer function in the digital control systems, are conducted in this article. Numerical integrators in such systems implement as an integral parts of the digital regulators. Objective. The goal of this research is an analysis of the behavior of explicit numerical integrators of different orders, which areused to discretize continuous systems, in order to study their impact on the properties of the synthesized digital system. Method. Numerical methods of integration are considered as digital filters, the behavior of which is studied by the frequency characteristics method. To do this, the z-transform apparatus was used. Integrators’ discrete transfer functions were found for frequency analysis using the Control Systems Toolbox package of the mathematical application MATLAB. For further analysis, two closed feedback test structures were used: with integrators in the forward channel and in the feedback loop. Both variants of structures were studied by the frequency characteristics method for sampling using numerical integrators of 1st–6th orders. Results. The inefficiency of using high-order numerical integrators for continuous systems’ discretization is shown. Given the behavior of the frequency characteristics of test systems, the most rational is the use of low-order integrators, namely – the first and second orders. Establishing the cause of this phenomenon requires additional research, in particular, to identify the possible impact of additional zeros and poles of discrete transfer functions of the numerical integrators. Conclusions. The use of low-order integrators, namely the first and second orders, is the most rational for sampling of digital control systems and the inefficiency of using high-order numerical integrators to sample continuous systems is proven.

Data-based stochastic modeling reveals sources of activity bursts in single-cell TGF-β signaling.

Cells sense their surrounding by employing intracellular signaling pathways that transmit hormonal signals from the cell membrane to the nucleus. TGF-β/SMAD signaling encodes various cell fates, controls tissue homeostasis and is deregulated in diseases such as cancer. The pathway shows strong heterogeneity at the single-cell level, but quantitative insights into mechanisms underlying fluctuations at various time scales are still missing, partly due to inefficiency in the calibration of stochastic models that mechanistically describe signaling processes. In this work we analyze single-cell TGF-β/SMAD signaling and show that it exhibits temporal stochastic bursts which are dose-dependent and whose number and magnitude correlate with cell migration. We propose a stochastic modeling approach to mechanistically describe these pathway fluctuations with high computational efficiency. Employing high-order numerical integration and fitting to burst statistics we enable efficient quantitative parameter estimation and discriminate models that assume noise in different reactions at the receptor level. This modeling approach suggests that stochasticity in the internalization of TGF-β receptors into endosomes plays a key role in the observed temporal bursting. Further, the model predicts the single-cell dynamics of TGF-β/SMAD signaling in untested conditions, e.g., successfully reflects memory effects of signaling noise and cellular sensitivity towards repeated stimulation. Taken together, our computational framework based on burst analysis, noise modeling and path computation scheme is a suitable tool for the data-based modeling of complex signaling pathways, capable of identifying the source of temporal noise.

A robust fully mixed finite element model for flow and transport in unsaturated fractured porous media

A fully mixed finite element (MFE) model is developed for nonlinear flow and transport in unsaturated fractured porous media with matrix-fracture and fracture-fracture fluid and mass exchanges. The model is based on the discrete fracture matrix (DFM) approach and assumes cross-flow equilibrium in the fractures. The MFE method is employed for the spatial discretization of both flow and transport on the 2D-matrix elements as well as on the 1D-fracture elements. An upwind scheme is employed to avoid unphysical oscillations in the case of advection dominant transport. The temporal discretization is performed using high-order time integration methods and efficient automatic time-stepping schemes via the MOL.Two test problems dealing with flow and mass transport in saturated and unsaturated fractured porous media are simulated to show the validity of the new model by comparison against (i) a 1D-2D Comsol finite element model and (ii) a 2D-2D Discontinuous Galerkin (DG) model where both fractures and matrix continua are discretized with small 2D mesh elements. The robustness and efficiency of the developed 1D-2D MFE model are then investigated for a challenging problem dealing with infiltration of contaminated water into an initially dry soil involving a fracture network.The new model yields stable results for advection-dominated and advection-dispersion transport configurations. Further, the results of the 1D-2D MFE model are in very good agreement with those of the 2D-2D DG model for both configurations. The simulation of infiltration of contaminated water into a dry fractured soil shows that the 1D-2D MFE model is within 15 times more efficient than the 2D-2D DG model, which confirms the high benefit of using robust and efficient DFM models for the simulation of flow and transport in fractured porous media.

Stable timestep formulas for high-order advection–diffusion and Navier–Stokes solvers

We derive three-dimensional stable timestep formulas for high-order explicit time integration of the advection–diffusion equation. The proposed formulas cover explicit first, second, third, and fourth-order Runge–Kutta integrators in time as well as upwind, central, second-order, high-order upwind (k-schemes), and flux-limiters for the advection term along with central discretization on the diffusion term. The stability rules are then compared to numerical computations of stable timesteps of the advection–diffusion equation and show excellent agreement. Finally, the formulas are implemented in a finite volume Navier–Stokes code and used to guide stable timestep selection to assess their suitability for Navier–Stokes solvers. All cases that were tested showed stable non-diverging solutions and significant improvement over common CFL-type stability specification. Despite its basis in the linear advection–diffusion equation, this work could be useful in providing a starting point for stable timestep selection for nonlinear evolutionary-type PDEs such as the Navier–Stokes equations.

Hamiltonian preserving nonlinear optics

In this paper we present a method of constructing a discrete nonlinear accelerator lattice that has an approximate integral of motion that is given upfront. The integral under consideration is a Hamiltonian in normalized (canonical) coordinates that is preserved by a lattice with a given accuracy. We establish a connection between the integrator of a Hamiltonian in normalized coordinates and a real lens arrangement. We consider accelerator as an analog computer and apply known algorithms of high-order symplectic integrators, to produce several nonlinear lattices. Following this concept we introduce new lattice design based on Ruth and Yoshida integrator that could be experimentally tested at existing accelerator facilities and Paul traps.

SIRIUS Project – IV. The formation history of the Orion Nebula Cluster driven by clump mergers

ABSTRACT The Orion Nebula Cluster (ONC) is an excellent example for understanding the formation of star clusters. Recent studies have shown that ONC has three distinct age populations and anisotropy in velocity dispersions, which are key characteristics for understanding the formation history of the ONC. In this study, we perform a smoothed-particle hydrodynamics/N-body simulation of star cluster formation from a turbulent molecular cloud. In this simulation, stellar orbits are integrated using a high-order integrator without gravitational softening; therefore, we can follow the collisional evolution of star clusters. We find that hierarchical formation causes episodic star formation that is observed in the ONC. In our simulation, star clusters evolve due to mergers of subclumps. The mergers bring cold gas with the clumps into the forming cluster. This enhances the star formation in the cluster centre. The dense cold gas in the cluster centre continues to form stars until the latest time. This explains the compact distribution of the youngest stars observed in the ONC. Subclump mergers also contribute to the anisotropy in the velocity dispersions and the formation of runaway stars. However, the anisotropy disappears within 0.5 Myr. The virial ratio of the cluster also increases after a merger due to the runaways. These results suggest that the ONC recently experienced a clump merger. We predict that most runaways originated from the ONC have already been found, but walkaways have not.

Energy-Preserving Continuous-Stage Exponential Runge--Kutta Integrators for Efficiently Solving Hamiltonian Systems

As one of the most important properties, energy preservation is a natural requirement for numerical integrators of Hamiltonian systems. Considering the limited second-order accuracy of the existing exponential average vector filed (or discrete gradient) method, this paper is devoted to developing high-order energy-preserving exponential integrators for general semilinear Hamiltonian systems. To this end, we first formulate and analyze the continuous-stage exponential Runge--Kutta (ERK) method. After deriving the order conditions, energy-preserving conditions, and symmetric conditions for the continuous-stage ERK method, we construct a novel class of fourth-order symmetric energy-preserving continuous-stage ERK methods with a free parameter based on the established conditions and present a convergence theorem for the fixed-point iteration associated with the continuous-stage ERK method. Furthermore, we extend the proposed continuous-stage ERK method to noncanonical Hamiltonian systems and discuss implementation issues in detail. Finally, numerical results from the FPU problem, the charged-particle dynamics, and the Klein--Gordon equation demonstrate the high efficiency, good energy-preserving behavior, and applicability of large time stepsizes of the proposed fourth-order energy-preserving continuous-stage ERK method.

Categorizing Visual Information in Subpopulations of Honeybee Mushroom Body Output Neurons.

Multisensory integration plays a central role in perception, as all behaviors usually require the input of different sensory signals. For instance, for a foraging honeybee the association of a food source includes the combination of olfactory and visual cues to be categorized as a flower. Moreover, homing after successful foraging using celestial cues and the panoramic scenery may be dominated by visual cues. Hence, dependent on the context, one modality might be leading and influence the processing of other modalities. To unravel the complex neural mechanisms behind this process we studied honeybee mushroom body output neurons (MBON). MBONs represent the first processing level after olfactory-visual convergence in the honeybee brain. This was physiologically confirmed in our previous study by characterizing a subpopulation of multisensory MBONs. These neurons categorize incoming sensory inputs into olfactory, visual, and olfactory-visual information. However, in addition to multisensory units a prominent population of MBONs was sensitive to visual cues only. Therefore, we asked which visual features might be represented at this high-order integration level. Using extracellular, multi-unit recordings in combination with visual and olfactory stimulation, we separated MBONs with multisensory responses from purely visually driven MBONs. Further analysis revealed, for the first time, that visually driven MBONs of both groups encode detailed aspects within this individual modality, such as light intensity and light identity. Moreover, we show that these features are separated by different MBON subpopulations, for example by extracting information about brightness and wavelength. Most interestingly, the latter MBON population was tuned to separate UV-light from other light stimuli, which were only poorly differentiated from each other. A third MBON subpopulation was neither tuned to brightness nor to wavelength and encoded the general presence of light. Taken together, our results support the view that the mushroom body, a high-order sensory integration, learning and memory center in the insect brain, categorizes sensory information by separating different behaviorally relevant aspects of the multisensory scenery and that these categories are channeled into distinct MBON subpopulations.

Contact-Implicit Direct Collocation With a Discontinuous Velocity State

Collisions between rigid bodies are impulsive events, meaning they create near-instantaneous, finite discontinuities in the velocity state of the dynamic system. If velocity is assumed to be continuous, that requirement may come into conflict with the desire for bodies to not interpenetrate, or with Coulomb friction (creating a condition known as the <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">Painlevé paradox).</i> This letter presents a contact-implicit framework for trajectory optimization using direct collocation that combines an unplanned impulsive contact model with an implicit, high-order numerical integration scheme. We demonstrate that this framework can resolve problems that other implicit formulations cannot, namely elastic or partially elastic collisions, and impacts without collision – the established resolution of the Painlevé paradox. We then evaluate its applicability to legged robotics problems and test its ability to discover the intricate sequence of varied collision types necessary to execute a skateboarding trick.

Distributed Prescribed-Time Consensus Observer for High-Order Integrator Multi-Agent Systems on Directed Graphs

This brief deals with the distributed consensus observer design problem for high-order integrator multi-agent systems on directed graphs, which intends to estimate the leader states accurately in a prescribed time interval. A new kind of distributed prescribed-time observers (DPTO) on directed graphs is first formulated for the followers, which is implemented in a cascading manner. Then, the prescribed-time zero-error estimation performance of the above DPTO is guaranteed for both time-invariant and time-varying directed interaction topologies, based on strictly Lyapunov stability analysis and mathematical induction method. Finally, the practicability and validity of this new distributed observer are illustrated via a numerical simulation example.

Application of Mayer’s activity expansions to the Ising problem

Recent successes in approximate evaluations of high-order cluster integrals have greatly advanced the quantitative applicability of Mayer’s cluster expansion for various statistical models of matter and even real substances. In this paper, the thermodynamics of the Ising model is represented in a general form of the virial series in powers of activity, which is not restricted by any special geometry or the simplification of nearest-neighbor interactions. In order to make the presented theoretical expressions applicable in quantitative studies, a number of techniques are considered, which approximate the unlimited (almost infinite) set of reducible cluster integrals, bn, based on only a few number of irreducible integrals, βk. Namely, such short βk sets are defined here for a number of magnetic models of various geometry, dimensionality, as well as interaction ranges (not limited by the nearest coordination sphere), and the calculations performed for these models in subcritical regimes (below the Curie point) indicate the adequacy of the proposed approach in general and good agreement with the known exact solution in a particular case of the two-dimensional square model with nearest-neighbor interactions. In addition, the quantitative similarity of subcritical magnetization curves is observed for all the considered models at the same values of a certain reduced temperature and the nature of this similarity may deserve further attentive studying.