Search

As the number of processors on supercomputers has increased dramatically, there is a growing interest in developing scalable algorithms with a high degree of parallelism for extreme-scale reservoir simulation. However, traditional simulators and algorithms for such nonlinear problems are usually based on the family of time-marching methods, where parallelization is restricted to the spatial dimension only. In this paper, we propose a family of parallel-in-time (PinT) reservoir algorithms for solving multiphase flow problems in porous media to fully exploit the parallelism of supercomputers, in which the design of efficient nonlinear and linear space-time algorithms plays an essential role. More precisely, we first introduce a space-time mixed finite element method for the fully-implicit discretization. The nonlinear system arising at each parallel-in-time step is solved by a variant of bound-preserving Newton methods with a nonlinear elimination preconditioner, where the corresponding linear system is solved by the space-time restricted additive Schwarz (stRAS) method. Since the straightforward extension of one-level stRAS to multilevel does not work because of the pollution effects from the temporal direction, we accordingly present a nonstandard V-cycle multilevel stRAS method with a pollution-removing strategy to accelerate the convergence and improve the robustness. Numerical experiments are presented to demonstrate that the aforementioned parallel-in-time algorithm can not only achieve a high degree of parallelism in accurately resolving reservoir transport features in a heterogeneous medium, but also successfully circumvent the convergence and scalability issues associated with the constraint of large time steps and the boundedness requirement of the solution. • Space-time mixed finite element scheme for parallel in-time discretization. • Active-set reduced-space method with nonlinear elimination for bound-preserving solution. • Space-time restricted additive Schwarz algorithm for linear preconditioning. • Nonlinear parallel-in-time simulations with up to 13,872 processor cores. [ABSTRACT FROM AUTHOR]

Machine learning methods are widely used in the natural sciences to model and predict physical systems from observation data. Yet, they are often used as poorly understood "black boxes," disregarding existing mathematical structure and invariants of the problem. Recently, the proposal of Hamiltonian Neural Networks (HNNs) took a first step towards a unified "gray box" approach, using physical insight to improve performance for Hamiltonian systems. In this paper, we explore a significantly improved training method for HNNs, exploiting the symplectic structure of Hamiltonian systems with a different loss function. This frees the loss from an artificial lower bound. We mathematically guarantee the existence of an exact Hamiltonian function which the HNN can learn. This allows us to prove and numerically analyze the errors made by HNNs which, in turn, renders them fully explainable. Finally, we present a novel post-training correction to obtain the true Hamiltonian only from discretized observation data, up to an arbitrary order. • The loss of Hamiltonian Neural Networks suffers from an artificial lower bound. • Using a symplectic integrator during training significantly improves performance. • Rigorous mathematical analysis yields an exact handle on the discretization error. • The learnt Hamiltonian can be corrected post-training, to any arbitrary order. [ABSTRACT FROM AUTHOR]

The optimal control of thermally convective flows is usually modeled by an optimization problem with constraints of Boussinesq equations that consist of the Navier-Stokes equation and an advection-diffusion equation. This optimal control problem is challenging from both theoretical analysis and algorithmic design perspectives. For example, the nonlinearity and coupling of fluid flows and energy transports prevent direct applications of gradient type algorithms in practice. In this paper, we propose an efficient numerical method to solve this problem based on the operator splitting and optimization techniques. In particular, we employ the Marchuk-Yanenko method leveraged by the L 2 -projection for the time discretization of the Boussinesq equations so that the Boussinesq equations are decomposed into some easier linear equations without any difficulty in deriving the corresponding adjoint system. Consequently, at each iteration, four easy linear advection-diffusion equations and two degenerated Stokes equations at each time step are needed to be solved for computing a gradient. Then, we apply the Bercovier-Pironneau finite element method for space discretization, and design a BFGS type algorithm for solving the fully discretized optimal control problem. We look into the structure of the problem, and design a meticulous strategy to seek step sizes for the BFGS efficiently. Efficiency of the numerical approach is promisingly validated by the results of some preliminary numerical experiments. [ABSTRACT FROM AUTHOR]

In this paper, a fully conservative sharp-interface method is developed for compressible multiphase flows with phase change. To enforce conservation, the coupling between two different phases is achieved via the interfacial fluxes, which are determined by solving a general Riemann problem with phase change. In addition, thermodynamical consistency is guaranteed by taking account of latent heat in the phase-change Riemann problem. This is efficiently and robustly solved by a novel approximate Riemann solver, where the complex eight-dimensional root-finding procedure required in the exact Riemann solver is simplified to solely iterating the mass flux. Compared to the existing methods, the present approximate Riemann solution consists of four waves, at which the jump conditions are imposed strictly to enforce the conservation of mass, momentum, and energy. Furthermore, for the evaluation of the mass flux, through numerical simulations and theoretical analysis, we have demonstrated that the adjacent states of a phase interface should be used to ensure numerical consistency, which, to our knowledge, has not been reported in the literature. Various numerical examples, including the phase-change Riemann problem, the oscillating droplet problem, and the shock-droplet interaction problem, are simulated to validate the present method. The numerical results are in good agreement with the exact solutions and the previous simulations. It is shown that, with the present method, the phase change effects in compressible multiphase flows are well resolved with conservation strictly preserved. • A fully conservative sharp-interface method is developed for compressible multiphase flows with phase change. • A novel approximate Riemann solver is developed to solve the two-phase phase-change Riemann problem. • The present Riemann solver employs a 4-wave model with all the jump conditions imposed strictly, which is more physical. • We ensure numerical consistency by using the adjacent states of the phase interface to compute the mass flux. [ABSTRACT FROM AUTHOR]

We assess the ability of three different approaches based on high-order discontinuous Galerkin methods to simulate under-resolved turbulent flows. The capabilities of the mass conserving mixed stress method as structure resolving large eddy simulation solver are examined. A comparison of a variational multiscale model to no-model or an implicit model approach is presented via numerical results. In addition, we present a novel approach for turbulent modeling in wall-bounded flows. This new technique provides a more accurate representation of the actual subgrid scales in the near wall region and gives promising results for highly under-resolved flow problems. In this paper, the turbulent channel flow and periodic hill flow problem are considered as benchmarks for our simulations. • Low dissipative high-order DG-FEM solver for implicit large eddy simulation. • Novel VMS like upwind stabilization for incompressible wall-bounded turbulent flows. • Comparison and verification using different high Reynolds number test cases. [ABSTRACT FROM AUTHOR]

In this paper, we are interested in constructing a scheme solving compressible Navier–Stokes equations, with desired properties including high order spatial accuracy, conservation, and positivity-preserving of density and internal energy under a standard hyperbolic type CFL constraint on the time step size, e.g., Δ t = O (Δ x). Strang splitting is used to approximate convection and diffusion operators separately. For the convection part, i.e., the compressible Euler equation, the high order accurate postivity-preserving Runge–Kutta discontinuous Galerkin method can be used. For the diffusion part, the equation of internal energy instead of the total energy is considered, and a first order semi-implicit time discretization is used for the ease of achieving positivity. A suitable interior penalty discontinuous Galerkin method for the stress tensor can ensure the conservation of momentum and total energy for any high order polynomial basis. In particular, positivity can be proven with Δ t = O (Δ x) if the Laplacian operator of internal energy is approximated by the Q k spectral element method with k = 1 , 2 , 3. So the full scheme with Q k (k = 1 , 2 , 3) basis is conservative and positivity-preserving with Δ t = O (Δ x) , which is robust for demanding problems such as solutions with low density and low pressure induced by high-speed shock diffraction. Even though the full scheme is only first order accurate in time, numerical tests indicate that higher order polynomial basis produces much better numerical solutions, e.g., better resolution for capturing the roll-ups during shock reflection. • The hyperbolic subproblem is solved by a positivity-preserving RKDG method. • The parabolic subproblem is solved by classical schemes to achieve positivity. • A proper choice of the IPDG method achieves total conservation. • Numerically a spatially higher order scheme produces better solutions. [ABSTRACT FROM AUTHOR]

The simulation of chemical kinetics involving multiple scales constitutes a modeling challenge (from ordinary differential equations to Markov chains) and a computational challenge (multiple scales, large dynamical systems, time step restrictions). In this paper, we propose a new discrete stochastic simulation algorithm: the postprocessed second kind stabilized orthogonal τ -leap Runge–Kutta method (PSK- τ -ROCK). In the context of chemical kinetics, this method can be seen as a stabilization of Gillespie's explicit τ -leap algorithm combined with a postprocessor. The stabilized procedure allows to simulate problems with multiple scales (stiff), while the postprocessing procedure allows to approximate the invariant measure (e.g. mean and variance) of ergodic stochastic dynamical systems. We prove stability and accuracy of the PSK- τ -ROCK for a reference system. Numerical experiments illustrate the high reliability and efficiency of the scheme when compared to other τ -leap methods. • Optimal explicit stabilized scheme for chemical kinetics. • Inexpensive ergodic accurate postprocessor. • Explicit scheme for stiff chemical kinetics problems. • Stabilized tau-leap method. [ABSTRACT FROM AUTHOR]

"Generalized Hydrodynamics" (GHD) stands for a model that describes one-dimensional integrable systems in quantum physics, such as ultra-cold atoms or spin chains. Mathematically, GHD corresponds to nonlinear equations of kinetic type, where the main unknown, a statistical distribution function f (t , z , θ) , lives in a phase space which is constituted by a one-dimensional position variable z , and a one-dimensional "kinetic" variable θ , actually a wave-vector, called "rapidity". Two key features of GHD equations are first a non-local and nonlinear coupling in the advection term, and second an infinite set of conserved quantities, which prevent the system from thermalizing. To go beyond this, we consider the dissipative GHD equations, which are obtained by supplementing the right-hand side of the GHD equations with a non-local and nonlinear diffusion operator or a Boltzmann-type collision integral. In this paper, we deal with new high-order numerical methods to efficiently solve these kinetic equations. In particular, we devise novel backward semi-Lagrangian methods for solving the advective part (the so-called Vlasov equation) by using a high-order time-Taylor series expansion for the advection fields, whose successive time derivatives are obtained by a recursive procedure. This high-order temporal approximation of the advection fields is used to design new implicit/explicit Runge–Kutta semi-Lagrangian methods, which are compared to Adams–Moulton semi-Lagrangian schemes. For solving the source terms, constituted by the diffusion and collision operators, we use and compare different numerical methods of the literature. • Numerical models for quantum experiments. • New models for quasi 1-d quantum systems. • Efficient numerical methods for such models. • Computer simulations for Bose Einstein Condensates. [ABSTRACT FROM AUTHOR]

Diffusion in heterogeneous media partitioned by semi-permeable interfaces has a wide range of applications in the physical and life sciences, including gas permeation in soils, diffusion magnetic resonance imaging (dMRI), drug delivery, thermal conduction in composite media, synaptic receptor trafficking, and intercellular gap junctions. At the single particle level, diffusion across a semi-permeable interface can be formulated in terms of so-called snapping out Brownian motion (SNOBM). The latter sews together successive rounds of reflected BM, each of which is restricted to one side of the interface. Each round of reflected BM is killed when the local time at the interface exceeds an independent, exponentially distributed random variable. (The local time specifies the amount of time a reflected Brownian particle spends in a neighborhood of the interface.) The particle then immediately resumes reflected BM on the same side or the other side of the interface according to a stochastic switch, and the process is iterated. In this paper, we develop a Monte Carlo algorithm for simulating a two-dimensional version of SNOBM, which is used to solve a first passage time (FPT) problem for diffusion in a domain with semi-permeable partially absorbing traps. Our method combines a walk-on-spheres (WOS) method with an efficient algorithm for computing the boundary local time that uses a Skorokhod integral representation of the latter. We validate our algorithm by comparing the Monte Carlo estimates of the mean first passage time (MFPT) to the exact solution for a single circular trap, and show that our simulations are consistent with asymptotic results obtained for the 2D narrow capture problem involving multiple small circular targets. We also discuss extensions to higher dimensions. • Developed a Monte Carlo algorithm for simulating diffusion through semi-permeable membranes. • Uses the probabilistic framework of snapping out Brownian motion. • Combines walk-on-spheres with an efficient method for computing the boundary local time. • Validated by comparing numerical results with analytic solutions of first passage time problems. [ABSTRACT FROM AUTHOR]

When solving compressible multi-material flow problems, an unresolved challenge is the computation of advective fluxes across material interfaces that separate drastically different thermodynamic states and relations. A popular idea in this regard is to locally construct bimaterial Riemann problems, and to apply their exact solutions in flux computation. For general equations of state, however, finding the exact solution of a Riemann problem is expensive as it requires nested loops. Multiplied by the large number of Riemann problems constructed during a simulation, the computational cost often becomes prohibitive. The work presented in this paper aims to accelerate the solution of bimaterial Riemann problems without introducing approximations or offline precomputation tasks. The basic idea is to exploit some special properties of the Riemann problem equations, and to recycle previous solutions as much as possible. Following this idea, four acceleration methods are developed, including (1) a change of integration variable through rarefaction fans, (2) storing and reusing integration trajectory data, (3) step size adaptation, and (4) constructing an R-tree on the fly to generate initial guesses. The performance of these acceleration methods is assessed using four example problems in underwater explosion, laser-induced cavitation, and hypervelocity impact. These problems exhibit strong shock waves, large interface deformation, contact of multiple (>2) interfaces, and interaction between gases and condensed matters. For all the problems, the acceleration methods are able to significantly reduce the computational cost without affecting solver robustness or solution accuracy. In different cases, the solution of bimaterial Riemann problems is accelerated by 37 to 87 times. As a result, the total cost of advective flux computation, which includes the exact Riemann problem solution at material interfaces and the numerical flux calculation over the entire computational domain, is accelerated by 18 to 81 times. • Bimaterial Riemann problems can be constructed to compute interfacial fluxes. • Computational cost is often prohibitive, if complex equations of state are involved. • An efficient exact Riemann solver is proposed, with four algorithmic improvements. • Significant cost reduction (37-87X) is achieved for four challenging examples. • Does not introduce approximations or require offline precomputation. [ABSTRACT FROM AUTHOR]

A fast and spectrally accurate Ewald summation method for the evaluation of stokeslet, stresslet and rotlet potentials of three-dimensional Stokes flow is presented. This work extends the previously developed Spectral Ewald method for Stokes flow to periodic boundary conditions in any number (three, two, one, or none) of the spatial directions, in a unified framework. The periodic potential is split into a short-range and a long-range part, where the latter is treated in Fourier space using the fast Fourier transform. A crucial component of the method is the modified kernels used to treat singular integration. We derive new modified kernels, and new improved truncation error estimates for the stokeslet and stresslet. An automated procedure for selecting parameters based on a given error tolerance is designed and tested. Analytical formulas for validation in the doubly and singly periodic cases are presented. We show that the computational time of the method scales like O (N log ⁡ N) for N sources and targets, and investigate how the time depends on the error tolerance and window function, i.e. the function used to smoothly spread irregular point data to a uniform grid. The method is fastest in the fully periodic case, while the run time in the free-space case is around three times as large. Furthermore, the highest efficiency is reached when applying the method to a uniform source distribution in a primary cell with low aspect ratio. The work presented in this paper enables efficient and accurate simulations of three-dimensional Stokes flow with arbitrary periodicity using e.g. boundary integral and potential methods. • A unified framework for fast summation of Stokes potentials in arbitrary periodicity. • Improved error estimates permit accurate and automatic selection of parameters. • New modified kernels ensure the truncation error to be independent of periodicity. • New analytical formulas permit validation in doubly and singly periodic cases. • Our method scales like O (N log ⁡ N) for N sources and target points. [ABSTRACT FROM AUTHOR]

In this paper, the new fifth-order, seventh-order, and ninth-order finite difference and finite volume multi-resolution weighted essentially non-oscillatory (MR-WENO) schemes with adaptive linear weights are presented for hyperbolic conservation laws on structured meshes. They are termed as high-order finite difference and finite volume ALW-WENO schemes. These ALW-WENO schemes only apply one small stencil and one large stencil in reconstruction processes, which could achieve the desired accuracy in the region of smoothness and non-oscillatory properties in the region of containing strong shocks. The linear weights that sum to one can be automatically adjusted to any positive numbers. This is the first time that arbitrary high-order finite difference and finite volume WENO schemes are designed by using only two unequal-sized central spatial stencils. The structure of these novel WENO schemes is simple, so it is easier for obtaining high-order accuracy and solving multi-dimensional problems in large scale engineering applications. Compared to traditional MR-WENO schemes with same order, the computational efficiency can be further improved. Some benchmark tests indicate that these new ALW-WENO schemes have good robustness and performance. • We propose a new class of MR-WENO schemes with adaptive linear weights. • The linear weights that sum to one can be automatically adjusted to any positive numbers. • We only use two unequal-sized central spatial stencils, which can achieve any high-order of accuracy in smooth regions. • The design process of these new WENO schemes is simpler and has higher computational efficiencies than before. • These new WENO schemes can be used for large scale engineering applications. [ABSTRACT FROM AUTHOR]

In this paper, we present a detailed numerical scheme for a single-phase compressible flow without diffusion of a multi-component mixture in porous media with the higher-order approximation in both space and time. The mathematical model consists of Darcy velocity, transport equations for each component of a mixture, pressure equation and associated relations for physical quantities such as viscosity or equation of state. The discrete problem is obtained using a combination of the discontinuous Galerkin method for the transport equations and the mixed-hybrid finite element method for the Darcy velocity and the pressure equation. In both methods the higher-order approximation is used. The resulting nonlinear problem for concentrations is solved with the fully mass-conservative iterative IMPEC method. Experimental order of convergence analysis (EOC) and some numerical experiments of a 2D flow are carried out. • A new numerical scheme is presented for simulation of single-phase flow in a porous medium. • The second-order approximation is attained in both space and time. • The second-order approximation in both time and space is verified using the experimental order of convergence. • Several computational examples illustrate advantageous features of our approach. [ABSTRACT FROM AUTHOR]

In this paper, we propose an arbitrary Lagrangian Eulerian (ALE) meshfree method that utilizes the radial basis function-finite difference (RBF-FD) method for solving diffusion-reaction equations on evolving surfaces. The surface evolution law is determined by a forced mean curvature flow (FMCF). We develop a parametric RBF-FD method to solving the FMCF, which has not been done previously in the framework of meshfree methods. The advantage of this parametric method is that it utilizes a suitable tangential velocity in the equation to maintain a fairly uniform distribution of nodes during the surface evolution. The tangential velocity also allows us to construct an ALE RBF-FD scheme for the surface PDEs, which performs better than the purely Lagrangian method. Furthermore, we discuss how to deal with the boundary conditions for both the FMCF and the PDEs on surfaces, and present a ghost node approach to efficiently discretize the boundary conditions. Some numerical experiments are shown, which not only demonstrate the effectiveness of the tangential velocity, but also illustrate two applications. [ABSTRACT FROM AUTHOR]

This paper develops and benchmarks an immersed peridynamics method to simulate the deformation, damage, and failure of hyperelastic materials within a fluid-structure interaction framework. The immersed peridynamics method describes an incompressible structure immersed in a viscous incompressible fluid. It expresses the momentum equation and incompressibility constraint in Eulerian form, and it describes the structural motion and resultant forces in Lagrangian form. Coupling between Eulerian and Lagrangian variables is achieved by integral transforms with Dirac delta function kernels, as in standard immersed boundary methods. The major difference between our approach and conventional immersed boundary methods is that we use peridynamics, instead of classical continuum mechanics, to determine the structural forces. We focus on non-ordinary state-based peridynamic material descriptions that allow us to use a constitutive correspondence framework that can leverage well-characterized nonlinear constitutive models of soft materials. The convergence and accuracy of our approach are compared to both conventional and immersed finite element methods using widely used benchmark problems of nonlinear incompressible elasticity. We demonstrate that the immersed peridynamics method yields comparable accuracy with similar numbers of structural degrees of freedom for several choices of the size of the peridynamic horizon. We also demonstrate that the method can generate grid-converged simulations of fluid-driven material damage growth, crack formation and propagation, and rupture under large deformations. [ABSTRACT FROM AUTHOR]

We introduce two new approximation methods for the numerical evaluation of the long-range component of the range-separated Coulomb potential and the approximation of the resulting high dimensional Two-Electron Integrals tensor (TEI) with long-range interactions arising in molecular simulations. The first method exploits the tensorized structure of the compressed two-electron integrals obtained through two-dimensional Chebyshev interpolation combined with Gaussian quadrature. The second method is based on the Fast Multipole Method (FMM). Numerical experiments for different medium size molecules on high quality basis sets outline the efficiency of the two methods. Detailed algorithmic is provided in this paper as well as numerical comparison of the introduced approaches. [ABSTRACT FROM AUTHOR]

Algebraic multigrid (AMG) is an effective iterative algorithm for solving large-scale linear systems. One challenge of constructing the AMG algorithm is the determination of the prolongation operator, which affects the convergence rate of AMG and is problem-dependent. In this paper, we propose a new Learning-based Local Weighted Least Squares (L-LWLS) method to construct the prolongation operator of AMG. Specifically, we construct the prolongation operator by solving the LWLS model with learned spatially-varying weights. We use the gradient descent algorithm to optimize the model with a learned initialization of the solution. Then the constructed prolongation operator is further corrected by a learned correction function to improve the convergence rate of AMG. We conduct experiments on solving graph Laplacian linear systems, diffusion partial differential equations, and Helmholtz equations. Experiments show that the proposed method can construct a better prolongation operator leading to a faster convergence rate than the compared methods, including the classical AMG, the smoothed aggregation AMG, the bootstrap AMG, and the learning-based AMG method. The results show that the proposed method can generalize well to different parameter distributions and problem sizes, i.e., the number of variables in the linear system. • We proposed a learning-based method L-LWLS to optimize the prolongation operator in AMG. • We propose to learn spatially-varying weights in LWLS and learn to correct the LWLS solution. • The proposed method accelerates the convergence rate of AMG for parameterized PDEs. • The proposed method can generalize to different parameter distributions and problem sizes. [ABSTRACT FROM AUTHOR]

This paper proposes high-order accurate well-balanced (WB) energy stable (ES) adaptive moving mesh finite difference schemes for the shallow water equations (SWEs) with non-flat bottom topography. To enable the construction of the ES schemes on moving meshes, a reformulation of the SWEs is introduced, with the bottom topography as an additional conservative variable that evolves in time. The corresponding energy inequality is derived based on a modified energy function, then the reformulated SWEs and energy inequality are transformed into curvilinear coordinates. A two-point energy conservative (EC) flux is constructed, and high-order EC schemes based on such a flux are proved to be WB that they preserve the lake at rest. Then high-order ES schemes are derived by adding suitable dissipation terms to the EC schemes, which are newly designed to maintain the WB and ES properties simultaneously. The adaptive moving mesh strategy is performed by iteratively solving the Euler-Lagrangian equations of a mesh adaptation functional. The fully-discrete schemes are obtained by using the explicit strong-stability preserving third-order Runge-Kutta method. Several numerical tests are conducted to validate the accuracy, WB and ES properties, shock-capturing ability, and high efficiency of the schemes. [ABSTRACT FROM AUTHOR]

This paper proposes a hierarchy of numerical fluxes for the compressible flow equations which are kinetic-energy and pressure equilibrium preserving and asymptotically entropy conservative, i.e., they are able to arbitrarily reduce the numerical error on entropy production due to the spatial discretization. The fluxes are based on the use of the harmonic mean for internal energy and only use algebraic operations, making them less computationally expensive than the entropy-conserving fluxes based on the logarithmic mean. The use of the geometric mean is also explored and identified to be well-suited to reduce errors on entropy evolution. Results of numerical tests confirmed the theoretical predictions and the entropy-conserving capabilities of a selection of schemes have been compared. [ABSTRACT FROM AUTHOR]

This paper studies three high-order structure-preserving finite volume weighted essentially non-oscillatory (WENO) methods, which are not only well balanced (WB) for a general known hydrostatic equilibrium state but also preserve the positivity of density and pressure, for the compressible Euler equations under gravitational fields. These methods are built on a simple local scaling positivity-preserving (PP) limiter and a modified WENO-ZQ reconstruction exactly preserving the cell average value and scaling invariance. The WB properties of these three methods are achieved based on suitable numerical fluxes and approximation to the gravitational source terms. Based on some convex decomposition techniques as well as several critical properties of the admissible states and numerical flux, we carry out rigorous positivity-preserving analyses for these three WB schemes. We rigorously prove that the three WB methods, coupled with the PP limiter and a strong-stability-preserving time discretization, are always PP under suitable Courant-Friedrichs-Lewy conditions. Extensive numerical examples are provided to confirm WB and PP properties of three methods. [ABSTRACT FROM AUTHOR]

This paper presents three Riemann solvers for the Wilkins model of hypoelasticity with the Jaumann derivative correction. The solvers are developed as extensions of the HLL and HLLC Riemann solvers for the Euler hydrodynamic model and differ in the number of waves considered (two-wave HLLEPJ Riemann solver, three-wave HLLCEPJ Riemann solver and five-wave HLLCEPJ Riemann solver). The discontinuous solution of the Wilkins system of equations is obtained using generalized Rankine-Hugoniot relations developed on the basis of the path-conservative DLM approach. A special choice of the path in the phase space allows to avoid any iterative processes in the solver construction. The proposed Riemann solvers are implicated in the first order Finite Volume numerical scheme for testing on the set of 1D, 2D and 3D problems of elastoplastic flows. • The non-iterative solvers differ in the number of waves considered - two, three and five waves are considered in the Riemann problem approximation. • The discontinuous solution of the Wilkins system of equations is obtained using generalized Rankine-Hugoniot relations developed on the basis of the path-conservative DLM approach. • The proposed solvers are tested by implication in the first order Finite Volume numerical scheme for numerical calculation of several test problems. [ABSTRACT FROM AUTHOR]

The Fully Implicit Method (FIM) is used widely in the reservoir simulation community, especially for the so-called black-oil model. The primary reason for resorting to FIM is its unconditional stability, but this comes at the high computational cost of solving coupled nonlinear systems of all the conservation equations simultaneously. In addition, having a scalable FIM simulator relies on specialized multi-stage linear solvers. In the last few years, the Sequential Fully Implicit (SFI) formulation has received a great deal of attention. SFI decomposes the fully coupled nonlinear problem into two nonlinear sub-problems: flow and transport. The flow problem is solved implicitly for pressure, and the transport problem is solved implicitly for the saturation(s). The splitting usually assumes that the total-velocity is fixed during the transport computations. Because the pressure has a near-elliptic behavior and the saturation(s) have near-hyperbolic behavior, the flow-transport splitting of SFI allows for the use of dedicated scalable nonlinear and linear solvers for each sub-problem. At convergence, SFI and FIM lead to the same numerical solution within the specified tolerances on the residuals. SFI, however, can suffer from convergence difficulties of the overall nonlinear iterations between flow and transport. The convergence difficulties of SFI are intimately related to the splitting, and they depend on the type and strength of coupling between the nonlinear flow and transport problems. In this paper, we develop indicators that quantify the nature and strength of coupling across the systems of equations and solution variables. The indicators quantify the coupling in both space and time. Our strategy employs terms in the FIM Jacobian matrix to quantifying the sensitivity of the mass conservation equations and the global volume-balance to changes in pressure and saturation. These 'coupling strength' indicators are used in a domain-decomposition strategy to define local sub-domains for which a fully coupled correction is performed instead of the usual fixed-point iteration update. We demonstrate the accuracy and utility of these indicators for several challenging problems, including buoyancy driven lock-exchange - with or without gas solubility, as well as, gas and water injection problems in heterogeneous reservoir models. • Local corrections in highly coupled regions help SFI with nonlinear convergence. • Coupling criteria help quantifying coupling strength and localizing coupled regions. • Mass balances sensitivity to pressure changes is capturing highly coupled subregions. [ABSTRACT FROM AUTHOR]

In this paper, a data-driven non-parametric approach is presented for forecasting the probability density evolution of stochastic dynamical systems. The method is based on stochastic Koopman operator and extended dynamic mode decomposition (EDMD). To approximate the finite-dimensional eigendecomposition of the stochastic Koopman operator, EDMD is applied to the training data set sampled from the stationary distribution of the underlying stochastic dynamical system. The family of the Koopman operators form a semigroup, which is generated by the infinitesimal generator of the stochastic dynamical system. A significant connection between the generator and Fokker-Planck operator is leveraged to construct an orthonormal basis of a weighted Hilbert space, enabling a spectral decomposition of the probability density function to be performed in this space. This approach is a data-driven method and used to predict the probability density evolution and real-time moment estimation. In the limit of the large number of snapshots and observables, the data-driven probability density approximation converges to the Galerkin projection of the semigroup solution of Fokker-Planck equation on a basis adapted to an invariant measure. This method shares similarities with diffusion forecasting but is capable of producing more accurate probability density forecasts than diffusion forecasting. A few numerical examples are presented to illustrate the performance of the proposed data-driven probability density forecast. • An orthonormal basis is constructed using the properties of infinitesimal generator of SDE and Fokker-Planck operator. • EDMD in weighted space is utilized to estimate the eigenpairs of stochastic Koopman operator. • The probability density function is approximated based on Koopman operator theory and spectral decomposition. • The proposed method is able to predict the statistical properties of nonlinear stochastic systems. [ABSTRACT FROM AUTHOR]

We present an open-source adaptive front-tracking solver for biological capsules in viscous flows. The membrane elastic and bending forces are solved on a Lagrangian triangulation using a linear Finite Element Method and a paraboloid fitting method. The fluid flow is solved on an octree adaptive grid using the open-source platform Basilisk. The Lagrangian and Eulerian grids communicate using an Immersed Boundary Method by means of Peskin-like regularized Dirac delta functions. We demonstrate the accuracy of our solver with extensive validations: in Stokes conditions against the Boundary Integral Method, and in the presence of inertia against similar (but not adaptive) front-tracking solvers. Excellent qualitative and quantitative agreements are shown. We then demonstrate the robustness of the present solver in challenging cases featuring extreme membrane deformations, very large computational domains and high volume fractions. Moreover, we illustrate the capability of the solver to simulate inertial capsule-laden flows in complex STL-defined geometries, opening the door for bioengineering applications featuring large three-dimensional channel structures. The source code and all the test cases presented in this paper are freely available. • Front Tracking/Immersed Boundary method to simulate a 3D flow laden with deformable capsules in a complex geometry. • Extension of the method to adaptive mesh refinement on Cartesian octree grids. • Implementation in the open-source software Basilisk. • Complex geometries are considered via STL files. • Comprehensive validations to verify the robustness and accuracy of the method. [ABSTRACT FROM AUTHOR]

An efficient and new methodology to deal with fluid structure interaction at high Reynolds number flows is presented in this article. It relies on the coupling of the lattice Boltzmann method and the immersed boundary method with a second order predictor corrector model for the structure. The effect of the lagrangian weight of the immersed boundary method is also analyzed in the context of fluid structure interaction. Both curved and moving boundaries are considered, and several methods to calculate the lagrangian weight are compared on relevant test-cases: the laminar flow around a cylinder at Reynolds number R e = 100 , the Poiseuille flow in a 2D channel and a 3D fluid-structure interaction test case: a deformable flapping flag immersed in a laminar flow. A convergence study in space and time is performed and the three following parameters are investigated: the number of lagrangian markers along the boundary, the value of the lagrangian weight and the shape of the discrete delta function used in the immersed boundary. Finally, the novelty of the paper is two-fold: a new expression of the lagrangian weight which is found to reduce the error by 20% in the case of fluid-structure interaction, and the coupling of a turbulence model with a second order predictor corrector model for improved stability and accuracy. This numerical methodology is found here to be accurate for challenging cases with high added mass effect and high Reynolds number flows. [ABSTRACT FROM AUTHOR]

In this paper, we present a novel kinetic scheme termed the Energy Conserving Semi-Lagrangian (ECSL) for the Vlasov-Ampère system. The novelty of the ECSL is that it retains the efficiency of the explicit scheme while simultaneously maintaining energy conservation and unconditional stability properties of the implicit scheme without relying on nonlinear iteration. The proposed ECSL method includes two main ingredients: the conservative Semi-Lagrangian (CSL) scheme and a novel field solver. The CSL scheme is utilized for the phase space discretization of the Vlasov equation, which enables the scheme to conserve mass exactly and remove the Courant-Friedrichs-Lewy restriction. The novel field solver is proposed by coupling the Ampère equation and the moments of the Vlasov equation in a semi-implicit way, allowing for an explicit and efficient calculation of the electric field. The combination of the novel field solver and CSL ensures that ECSL scheme conserves total energy and mass on the fully discrete level, regardless of spatial resolution and time step size. Moreover, the ECSL scheme still provides reliable solution even when its spatial and temporal resolution are insufficient to fully resolve the Debye length and plasma period, making it a promising tool for multiscale and lengthy simulations. Several numerical experiments are presented to demonstrate the accuracy, efficiency, and conservation properties of the proposed method. • We present a novel energy conserving semi-Lagrangian (ECSL) method for electrostatic plasma simulations. • ECSL retains the efficiency and simplicity of the explicit scheme while conserves energy exactly. • ECSL is unconditionally stable, freeing the user from the need to resolve the Debye length and plasma period. • The first grid-based kinetic method that attains all the above properties. • Higher orders in both space and time and easily extended to multidimensional simulations. [ABSTRACT FROM AUTHOR]

This paper proposes a novel computational modeling approach to investigate the fluid-structure interactions with moving contact lines. By embracing the generalized Onsager principle, a coupled hydrodynamics and phase field system is introduced that can describe the fluid-structure interactions with moving contact lines in a thermodynamically consistent manner. The resulting partial differential equation (PDE) model consists of the Navier-Stokes equation and a nonlinear Allen-Cahn type equation. Volume conservation is enforced through an additional penalty term. A fully-discrete structure-preserving numerical scheme is proposed by combining several techniques to solve this coupled PDE system effectively and accurately. For the temporal discretization, we utilize the supplementary variable method for preserving the thermodynamic structure and the projection approach for reducing the problem size. Then, we use the finite difference method on the staggered grid for spatial discretization. Furthermore, we have rigorously proved that the proposed numerical scheme based on the second-order backward difference formula respects the original energy stability, i.e., the scheme is energy stable. Additionally, with the aid of the supplementary variable method, the resultant scheme can be transformed into a constrained optimization problem, where the solutions of the supplementary variables are the arguments of the objective function that reaches the optimality. Then the augmented Lagrangian method is introduced in part to bring robustness and efficiency to solving such a constrained optimization problem. Finally, various numerical simulations verify the model's capability and demonstrate the scheme's effectiveness, accuracy, and stability. • We propose a novel model for fluid-structure interactions with moving contact lines. • Our PDE models embrace the generalized Onsager principle for thermodynamic consistency. • We introduce a supplementary variable to change the numerical scheme design into an optimization problem. • Contrary to other schemes, our numerical approach respects the original energy dissipation law. [ABSTRACT FROM AUTHOR]

In this paper, we introduce an approximation method that establishes certain order enhancements by leveraging radial basis functions (RBFs) in the numerical solution of conservation laws. The use of RBFs for interpolation and approximation is a well developed area of research. Of particular interest in this work is the development of high order finite volume (FV) weighted essentially non-oscillatory (WENO) methods, which utilize RBF approximations to obtain required data at cell interfaces. The aforementioned improvement in the order of accuracy is addressed through an analysis of the truncation error, resulting in expressions for the shape parameters appearing in the basis. This paper seeks to address the practical elements of the approach, including the evaluations of shape parameters as well as a hybrid implementation. To highlight the effectiveness of the non-polynomial basis in shock-capturing, the proposed methods are applied to systems of one-dimensional hyperbolic and weakly hyperbolic conservation laws and compared with several well-known WENO schemes in the literature. We also include a two-dimensional example for a scalar problem that demonstrates an extension to multiple dimensions. In the case of the non-smooth, weakly hyperbolic test problem, notable improvements are observed in predicting the location and height of the finite time blowup. The numerical results demonstrate that the proposed schemes attain notable improvements in accuracy, as indicated by the analysis of the reconstructions. A key contribution of this work is the development of robust third-order WENO method, which further demonstrates the effectiveness of the non-polynomial basis. • Order enhanced FV WENO methods which utilize radial basis function(RBF) approximation. • Order of accuracy is addressed through an analysis for the shape parameters in RBFs. • The implementation employs a hybrid strategy, aims to alleviate computational costs. • Applied to hyperbolic and weakly hyperbolic conservation laws with well-known schemes. [ABSTRACT FROM AUTHOR]

Simulating the dynamics of a nonequilibrium quantum many-body system by computing the two-time Green's function associated with such a system is computationally challenging. However, we are often interested in the time-diagonal of such a Green's function or time-dependent physical observables that are functions of one time. In this paper, we discuss the possibility of using dynamic mode decomposition (DMD), a data-driven model order reduction technique, to characterize one-time observables associated with the nonequilibrium dynamics using snapshots computed within a small time window. The DMD method allows us to efficiently predict long time dynamics from a limited number of trajectory samples. We demonstrate the effectiveness of DMD on a model two-band system. We show that, in the equilibrium limit, the DMD analysis yields results that are consistent with those produced from a linear response analysis. In the nonequilibrium case, the extrapolated dynamics produced by DMD is more accurate than a special Fourier extrapolation scheme presented in this paper. We point out a potential pitfall of the standard DMD method caused by insufficient spatial/momentum resolution of the discretization scheme. We show how this problem can be overcome by using a variant of the DMD method known as higher order DMD. • Perform accurate and efficient simulation of nonequilibrium many-body dynamics. • Use Dynamic mode decomposition (DMD) to approximate two-time Green's function. • Improve the accuracy of DMD by higher order DMD (HODMD). • Demonstrate the efficiency and accuracy of DMD and HODMD on a model two-band system. [ABSTRACT FROM AUTHOR]

With the increasingly accurate modeling and simulation demands and techniques, obtaining the knowledge of the multi-component phase equilibrium states holds the bottleneck challenge in the study of multi-phase fluid flow. In order to resolve the shortages in computational efficiency and stability in the conventional iterative flash calculation, a new phase equilibrium prediction method is proposed by solving the saddle point of the multi-phase system. In this paper, the HiSD (High-index saddle point dynamics) algorithm is used for the first time to calculate the saddle points on the energy landscape of the two-phase two-component NVT flash model based on the Peng-Robinson equation of state, and the up-down search algorithm of HiSD is applied to generate the solution landscape of the flash system. The Rosenbrock-Euler ETD (exponential time differencing) format is involved to reduce the interference of system rigidity to the calculation. It can be referred from the numerical analysis that there are at most the 1 s t -order saddle points in the energy landscape of the two-component two-phase NVT flash system, and all these saddle points are located on one straight line of the hyperplane, where the energy is equal everywhere. All these 1 s t -order saddle points can converge to the same or equivalent local minima, which indicates that the two-component two-phase flash system is a system with only one single solution with physical meanings. In addition, the saddle points also obey a linear relationship and the energy remains the same at different temperatures. Therefore, using the method proposed in this paper, the conventional two-step efforts of phase stability test and phase separation calculation can be simplified. The 1 s t -order saddle points of the system can be directly calculated, reducing the need for an initial guess. The local minima can be directly searched through the downward direction of the saddle point, which greatly reduces the calculation amount of phase equilibrium calculations. Furthermore, the minimum states at different temperatures can be calculated in batch by using one certain initial value, which significantly improves the adaptability and reliability for complex engineering problems with drastic temperature changes. • The HISD algorithm used for the first time in saddle point calculation in NVT flash. • The Rosenbrock-Euler ETD format used to reduce the interference of system rigidity. • At most the 1st-order saddle points in the energy landscape of the 2-component system. • There is only one solution in the 2-component 2-phase NVT flash problem. • All the 1st-order saddle points can converge to the equivalent local minima. [ABSTRACT FROM AUTHOR]

This paper presents the first parallel implementation of the novel "Interpolated Factored Green Function" (IFGF) method introduced recently for the accelerated evaluation of discrete integral operators arising in wave scattering and other areas (Bauinger and Bruno, Jour. Computat. Phys., 2021). On the basis of the hierarchical IFGF interpolation strategy, the proposed (hybrid MPI-OpenMP) parallel implementation results in efficient data communication, and it scales up to large numbers of cores—without any hard limitations on the number of cores efficiently employed. Moreover, on any given number of cores, the proposed parallel approach preserves the O (N log ⁡ N) computing cost inherent in the sequential version of the IFGF algorithm. Unlike other approaches, the IFGF method does not utilize the Fast Fourier Transform (FFT), and it is thus better suited for efficient parallelization in distributed-memory computer systems. In particular, the IFGF method relies on a "peer-to-peer" strategy wherein, at every level, field propagation is directly enacted via "exchanges" between "peer" polynomials of constant degree, without data accumulation in large-scale "telephone-central" mathematical constructs such as those used in the Fast Multipole Method (FMM) and pure FFT-based approaches. A variety of numerical results presented in this paper illustrate the character of the proposed parallel algorithm, in particular demonstrating scaling from 1 to all 1,680 cores available in the High Performance Computing cluster used, and for problems of up to 4,096 wavelengths in acoustic size. • New interpolation-based Green function acceleration method: IFGF. • Effective OpenMP-MPI parallelization strategy, small memory footprint. • Avoidance of parallelization bottlenecks, high strong parallel scaling efficiency. • Demonstration on geometries up to 4,096 wavelengths in size. • Acceleration of discrete operators containing over 2 billion degrees of freedom. [ABSTRACT FROM AUTHOR]

In this paper, we evaluate the performance of the multilevel Monte Carlo method (MLMC) for deterministic and uncertain hyperbolic systems, where randomness is introduced either in the modeling parameters or in the approximation algorithms. MLMC is a well known variance reduction method widely used to accelerate Monte Carlo (MC) sampling. However, we demonstrate in this paper that for hyperbolic systems, whether MLMC can achieve a real boost turns out to be delicate. The computational costs of MLMC and MC depend on the interplay among the accuracy (bias) and the computational cost of the numerical method for a single sample, as well as the variances of the sampled MLMC corrections or MC solutions. We characterize three regimes for the MLMC and MC performances using those parameters, and show that MLMC may not accelerate MC and can even have a higher cost when the variances of MC solutions and MLMC corrections are of the same order. Our studies are carried out by a few prototype hyperbolic systems: a linear scalar equation, the Euler and shallow water equations, and a linear relaxation model, the above statements are proved analytically in some cases, and demonstrated numerically for the cases of the stochastic hyperbolic equations driven by white noise parameters and Glimm's random choice method for deterministic hyperbolic equations. • Investigate multilevel Monte Carlo for deterministic and uncertain hyperbolic systems. • Identify three regimes for the MLMC/MC performances using bias, cost and variance. • MLMC may be outperformed by MC when their variances are of the same order. • Analytical proof for prototypical scalar advection equation with white noise parameters. • Numerical justification for stochastic hyperbolic equations and Glimm's method. [ABSTRACT FROM AUTHOR]

In this paper, we introduce a new reduced basis methodology for accelerating the computation of large parameterized systems of high-fidelity integral equations. Core to our methodology is the use of coarse-proxy models (i.e., lower resolution variants of the underlying high-fidelity equations) to identify important samples in the parameter space from which a high quality reduced basis is then constructed. Unlike the more traditional POD or greedy methods for reduced basis construction, our methodology has the benefit of being both easy to implement and embarrassingly parallel. We apply our methodology to the under-served area of integral equations, where the density of the underlying integral operators has traditionally made reduced basis methods difficult to apply. To handle this difficulty, we introduce an operator interpolation technique, based on random sub-sampling, that is aimed specifically at integral operators. To demonstrate the effectiveness of our techniques, we present two numerical case studies, based on the radiative transport equation and a boundary integral formation of the Laplace equation respectively, where our methodology provides a significant improvement in performance over the underlying high-fidelity models for a wide range of error tolerances. Moreover, we demonstrate that for these problems, as the coarse-proxy selection threshold is made more aggressive, the approximation error of our method decreases at an approximately linear rate. Finally, we provide a public repository of our source code with easy instructions for reproducing all results in this paper. • Presents novel model reduction method for many-query integral equation problems. • Uses column-pivoted QR and coarse proxy models to efficiently extract training samples for a reduced basis. • Sparse sampling of integral operators allows for efficient reconstruction of Galerkin projected operators. • Benchmarks on radiative transport and Laplace equation show a 30× speedup with only a 4% sacrifice to accuracy. [ABSTRACT FROM AUTHOR]

On the occasion of 60 years to a seminal paper by the Swedish mathematician Germund Dahlquist, this review paper starts by discussing the celebrated Dahlquist's barriers , which are theoretical limitations on the accuracy and stability properties of a broad and important class of time-stepping methods, i.e., Linear Multistep (LMS) methods, for the solution of initial value problems. Perhaps the two most dramatic Dahlquist barriers are the one that precludes an explicit LMS method from being unconditionally stable, and the one that precludes a high-order LMS method from being unconditionally stable. We discuss Dahlquist's barriers and also later barriers proved by other authors. We then discuss some time-stepping methods which seemingly break at least one of Dahlquist's barriers. Of course, the explanation of this "paradox" is that these are not LMS methods, so the barriers do not necessarily apply to them. We relate to two types of barrier breakers: those that break stability barriers (explicit unconditionally-stable methods), and those that break order barriers (high-order explicit or unconditionally-stable methods). We also review some current advanced time-integration techniques which significantly deviate from the LMS format. • The Dahlquist barriers are presented in their full beauty. • Other barriers for time-integration schemes are discussed. • Methods which "break" the stability barriers are discussed. • Methods which "break" the order barriers are discussed. • Recent advanced methods for time integration are reviewed. [ABSTRACT FROM AUTHOR]

In this paper, a family of stencil selection algorithms is presented for WENO schemes on unstructured meshes. The associated freedom of stencil selection for unstructured meshes, in the context of WENO schemes present a plethora of various stencil selection algorithms. The particular focus of this paper is to assess the performance of various stencil selection algorithm, investigate the parameters that dictate their robustness, accuracy and computational efficiency. Ultimately, efficient and robust stencils are pursued that can provide significant savings in computational performance, while retaining the non-oscillatory character of WENO schemes. This is achieved when making the stencil selection algorithms adaptive, based on the quality of the cells for unstructured meshes, that can in turn reduce the computational cost of WENO schemes. For assessing the performance of the developed algorithms well established test problems are employed. These include the least square approximation of polynomial functions, linear advection equation of smooth functions and solid body rotation test problem. Euler and Navier-Stokes equations test problems are also pursued such as the Shu-Osher test problem, the Double Mach Reflection, the supersonic Forward Facing step, the Kelvin-Helmholtz instability, the Taylor-Green Vortex, and the flow past a transonic circular cylinder. • Various central and directional stencils algorithms presented. • Uniqueness, Cover, Width and Compactness the most dominant parameters for directional stencils. • Algorithms suitability is dependent on the numerical framework. • Adaptive stencils reduce the computational cost of the WENO schemes. [ABSTRACT FROM AUTHOR]

In this paper, we are interested in propagating uncertainties through the linear Boltzmann equation. Such model is intensively used in neutronics, photonics, socio-economics, epidemiology. It is often solved thanks to Monte-Carlo (MC) schemes. MC codes are reliable and accurate but costly and, as a consequence, propagating uncertainties through them remains quite complicated. In order to alleviate the cost of an uncertainty propagation, reduced models are often built. In this paper, we focus on generalised Polynomial Chaos (gPC) reduced models, and especially on their resolution with an MC scheme: such strategy is commonly called Monte Carlo-generalised Polynomial Chaos (MC-gPC) in the literature [1–10]. It allows important computational gains on many applications: in a nutshell, the reasons for its success are spectral convergence [2] plus the fact that it is based on an MC resolution which can be implemented thanks to simple modifications of an existing MC code [1,9,8]. But MC-gPC also presents some weaknesses: it is sensitive to the curse of dimensionality [1,9] and is noisier than other strategies [10]. The aim of this paper is to present new MC schemes solving the same gPC based reduced model but attenuating the two previous drawbacks. They are based on multigroup-like resolution methods. The new MC schemes improve the run times of MC-gPC. The resolution scheme is intrusive: this means that modifications of an existing solver are necessary (even if people familiar with multigroup MC resolution will not be intimidated by them). The paper ends with a discussion about taking into account uncertainties at the early stages of the development of a simulation code together with some original and efficient hybrid intrusive/non-intrusive applications. • New solver of Intrusive generalised Polynomial Chaos (gPC). • More efficient than non-intrusive gPC and more efficient than former MC-gPC versions. • Monte Carlo (MC) resolution of the gPC based reduced model with a multigroup-like scheme. • Efficient uncertainty propagation with a combined MC/gPC Chaos resolution. • Efficient Hybrid non-intrusive/intrusive application. [ABSTRACT FROM AUTHOR]

This paper presents under-resolved simulations of compressible flows in complex structures like buildings. The mean pressure field is the primary focus as it is of high importance in many cases, such as the overpressure generated by the detonation of an explosive charge. In such high-risk situations, the mean pressure in the different building rooms has to be determined as quickly as possible, in order to help first responders evaluating the risk level of entering the building right after an explosion for instance. Conventional computational fluid dynamics methods are obviously able to deliver those requested pressure fields. However, they also require a long and tedious meshing process due to the presence of small openings in the computational domain (doors, windows or stairwells). In turn, these small openings need spatial and temporal resolution having significant consequences on computation time. But when dealing with pressing situations, a long pre-processing is not acceptable. The present paper proposes a two-step solution method. First, a very rough mesh is constructed which does not consider small openings. The mesh is constructed in an easy and particular way such that simplifications can be made when computing the fluid flux through these openings. Second, a specific Riemann solver dealing with geometric discontinuities is developed to estimate the fluid flow. The Riemann solver takes advantage of the simplifications provided by the meshing method, and accounts for both unchoked and choked flows through openings. The overall method then consists of both a specific meshing process that requires very little pre-processing, which considerably reduces the time needed to produce a complete simulation, and a Riemann solver. The proposed method is named MUZO in reference to its "MUlti-ZOne" flow solver. The MUZO solver is validated against resolved computations and provides accurate simulated mean pressure fields within a realistic building in a few seconds. • 3D simulations of compressible flow discharge in buildings are addressed through a very fast numerical framework. • Only 2D footprints without geometric details (doors, windows, etc.) are required to generate 3D meshes. • The previously omitted geometric restrictions are considered through a specific Riemann solver. [ABSTRACT FROM AUTHOR]

In this paper, we construct sign-preserving second-order IMplicit Pressure Explicit Concentration (IMPEC) time methods for generalized coupled non-Darcy flow and transport problems in petroleum engineering, extending the algorithm given in [10] which is only applicable to Darcy flows. We use interior penalty discontinuous Galerkin (IPDG) methods for spatial discretization, and develop bound-preserving technique to obtain physically relevant numerical approximations. The sign-preserving second-order IMPEC method is an important innovation. The method is based on the framework of the second-order strong-stability-preserving Runge-Kutta (SSP-RK2) method. The basic idea is to treat the pressure equation implicitly and the concentration equation explicitly so as to obtain a first-order time integration. Then we introduce a correction stage to compensate the accuracy, maintaining the physical bounds of the numerical cell averages. Unfortunately, the above algorithm is not applicable to non-Darcy problems. There are two main difficulties. Firstly, since the velocity equation is nonlinear and time-independent, all variables in the equation must be calculated at the same time level. However, the treatment in [10] will yield an extremely complicated algorithm and significantly large computational cost, and some iterations whose convergence may not be available if the solutions are not smooth. In our scheme, we linearize the velocity equation and the numerical solutions are not at the same time level, leading to first-order accurate solutions. Therefore, we adopt a completely different approach from [10] to derive the correction stage, keeping the physical bounds of the numerical solutions. Secondly, though with the new correction stage, it is still not easy to solve velocity equations in some stages. In this paper, we construct a direct solver for the velocity equation to save computational cost. Numerical experiments will be given to demonstrate that the improved sign-preserving second-order IMPEC scheme can reduce the computational cost significantly compared with explicit schemes if the diffusion coefficient D is small in the concentration equation. The proposed method also yields much larger CFL number compared with first-order IMPEC schemes. Moreover, the effectiveness of the bound-preserving technique will also be verified. • Construct second-order implicit pressure explicit concentration methods. • The algorithm works for both Darcy and non-Darcy flows. • The method saves more than 90% of the CPU time. • The algorithm preserves the physical bounds of the concentration. • Provide a fast solver for the Darcy-Forchheimer model without iterations. [ABSTRACT FROM AUTHOR]

In this paper, we present a general derivation of multicontinuum equations and discuss cell problems. We present constraint cell problem formulations in a representative volume element and oversampling techniques that allow reducing boundary effects. We discuss different choices of constraints for cell problems. We present numerical results that show how oversampling reduces boundary effects. Finally, we discuss the relation of the proposed methods to our previously developed methods, Nonlocal Multicontinuum Approaches. • We derive multicontinuum methods using a homogenization-like expansion and present constraint cell problem formulations. • Constraint cell problems allow using averages for different continua and give a flexibility to the framework. • We discuss appropriate local boundary conditions in representative volume elements and introduce oversampling. • The resulting multicontinuum equations show that local averages of the solution will differ among each other. • The average constraints, discussed in this paper, are easy to set and guarantee exponential decay. [ABSTRACT FROM AUTHOR]

In the framework of fusion energy research, divertor design and model calibration based on plasma edge codes currently rely either on manual iterative tuning of parameters, or alternatively on parameter scans. This, combined with the complex and computationally expensive nature of plasma edge codes, makes these procedures extremely cumbersome. Gradient-based optimization methods can significantly reduce this effort, but require an efficient strategy for sensitivity calculation. Algorithmic Differentiation (AD) offers an efficient and accurate solution, allowing semi-automatic sensitivity computation in complex, continuously developed codes. The adjoint AD mode is especially attractive, as its cost is independent of the number of input parameters. In this paper, adjoint AD sensitivity calculation is deployed for the first time in plasma edge codes, applying the TAPENADE AD tool to SOLPS-ITER. Adjoint AD results are verified to be machine precision accurate compared to tangent AD mode, and up to 10−9 compared to finite differences. Scalings of AD computational efforts prove the advantages of adjoint compared to tangent AD, while memory requirements rapidly increase for adjoint, showing the need for an improved checkpointing strategy. With the new tool, a wealth of information becomes accessible to the modeler. An adjoint sensitivity analysis for a COMPASS density scan identifies the input parameters with largest impact on the solution, and 2D sensitivity maps show their spatial dependence. • Algorithmic Differentiation allows adjoint sensitivity computation in complex codes. • Algorithmic Differentiation is applied to the plasma edge code SOLPS-ITER. • Adjoint sensitivities are obtained for the first time in plasma edge codes. • The computational cost of the sensitivities agrees with the expected performance. • A large amount of new sensitivity information is made available for the community. [ABSTRACT FROM AUTHOR]

Jump penalty stabilization techniques for under-resolved turbulence have been recently proposed for continuous and discontinuous high order Galerkin schemes [1–3]. The stabilization relies on the gradient or solution discontinuity at element interfaces to incorporate localised numerical diffusion in the numerical scheme. This diffusion acts as an implicit subgrid model and stabilizes under-resolved turbulent simulations. This paper investigates the effect of jump penalty stabilization methods (penalising gradient or solution) for stabilization and improvement of high-order discontinuous Galerkin schemes in turbulent regime. We analyze these schemes using an eigensolution analysis, a 1D non-linear Burgers equation (mimicking a turbulent cascade) and 3D turbulent Navier-Stokes simulations (Taylor-Green Vortex and Kelvin-Helmholtz Instability problems). We show that the two jump penalty stabilization techniques can stabilize under-resolved simulations thanks to the improved dispersion-dissipation characteristics (when compared to non-penalized schemes) and provide accurate results for turbulent flows. The numerical results indicate that the proposed jump penalty methods stabilize under-resolved simulations and improve the simulations, when compared to the original unpenalized scheme and to classic explicit subgrid models (Smagorinsky and Vreman). • Stabilisation of DG for under-resolved turbulence through penalising gradient or solution jumps at element interfaces. • Implementation for 1D advection-diffusion and 3D Navier-Stokes equations. • Eigensolution analysis for dispersion-dissipation behaviour for various polynomial orders and Riemann fluxes. • Validation on Burgers and Taylor-Green vortex problems and comparison with classic explicit subgrid models. [ABSTRACT FROM AUTHOR]

The opacity of FRTE depends on not only the material temperature but also the frequency, whose values may vary several orders of magnitude for different frequencies. The gray radiation diffusion and frequency-dependent diffusion equations are two simplified models that can approximate the solution to FRTE in the thick opacity regime. The frequency discretization for the two limit models highly affects the numerical accuracy. However, classical frequency discretization for FRTE considers only the absorbing coefficient. In this paper, we propose a new decomposed multi-group method for frequency discretization that is not only AP in both gray radiation diffusion and frequency-dependent diffusion limits, but also the frequency discretization of the limiting models can be tuned. Based on the decomposed multi-group method, a full AP scheme in frequency, time, and space is proposed. Several numerical examples are used to verify the performance of the proposed scheme. • The frequency discretization can be consistent with the required multi-group diffusion model. • Accuracy is guaranteed in both optically thin and thick regimes with large spatial meshes. • The required temporary step is independent of the light speed. • Nonlinearity appears only in the system with macroscopic quantities. • The microscopic quantities can be updated linearly. [ABSTRACT FROM AUTHOR]

Recently, there has been a growing interest in efficient numerical algorithms based on tensor networks and low-rank techniques to approximate high-dimensional functions and solutions to high-dimensional PDEs. In this paper, we propose a new tensor rank reduction method based on coordinate transformations that can greatly increase the efficiency of high-dimensional tensor approximation algorithms. The idea is simple: given a multivariate function, determine a coordinate transformation so that the function in the new coordinate system has smaller tensor rank. We restrict our analysis to linear coordinate transformations, which gives rise to a new class of functions that we refer to as tensor ridge functions. Leveraging Riemannian gradient descent on matrix manifolds we develop an algorithm that determines a quasi-optimal linear coordinate transformation for tensor rank reduction. The results we present for rank reduction via linear coordinate transformations open the possibility for generalizations to larger classes of nonlinear transformations. • Riemannian optimization of tensor ridge functions. • Tensor rank reduction via coordinate flows. • PDEs in curvilinear coordinate systems. • Rank-adaptive tensor methods. [ABSTRACT FROM AUTHOR]

Summation-by-parts (SBP) operators allow us to systematically develop energy-stable and high-order accurate numerical methods for time-dependent differential equations. Until recently, the main idea behind existing SBP operators was that polynomials can accurately approximate the solution, and SBP operators should thus be exact for them. However, polynomials do not provide the best approximation for some problems, with other approximation spaces being more appropriate. We recently addressed this issue and developed a theory for one-dimensional SBP operators based on general function spaces, coined function-space SBP (FSBP) operators. In this paper, we extend the theory of FSBP operators to multiple dimensions. We focus on their existence, connection to quadratures, construction, and mimetic properties. A more exhaustive numerical demonstration of multi-dimensional FSBP (MFSBP) operators and their application will be provided in future works. Similar to the one-dimensional case, we demonstrate that most of the established results for polynomial-based multi-dimensional SBP (MSBP) operators carry over to the more general class of MFSBP operators. Our findings imply that the concept of SBP operators can be applied to a significantly larger class of methods than is currently done. This can increase the accuracy of the numerical solutions and/or provide stability to the methods. • Multi-dimensional summation-by-parts operators for general (non-polynomial) function spaces. • Existence proof that connects them to quadratures. • Construction procedure for general multi-dimensional domains/elements/blocks. • Proof of discrete conservation and energy-boundedness. • Numerical comparison demonstrating that non-polynomial function spaces can be of advantage. [ABSTRACT FROM AUTHOR]

The discrete-equations method (DEM) [1] provides a universal approach to solve multi-phase-flow equations as it combines the solutions of pairwise Riemann problems. Although very robust, the original DEM with piecewise-constant volume fractions suffers from strong diffusion preventing accurate interface capturing. High-order interface reconstruction, however, introduces a restrictive time-step limit. This paper presents RDEMIC, a robust extension of DEM for accurate interface capturing on Cartesian meshes. By a modified partitioning of the Riemann solutions and a specific combination of fluxes and non-conservative terms, the time-step restriction is effectively prevented, which is critical for making the method practically applicable. Moreover, the accuracy of interface and shock-wave propagation is maintained. RDEMIC is not limited to two-phase flow but defined for an arbitrary number of phases. The method is combined with a THINC scheme [2] to reconstruct volume fractions. The reconstruction is enhanced by a positivity-preserving averaging procedure, which is consistent with the underlying multi-stage Runge–Kutta scheme of the flow solver. The resulting scheme consisting of RDEMIC and the positivity-preserving THINC reconstruction is very robust and captures the interface with high accuracy. We demonstrate its performance for various cases of shock-interface interactions, which show very good agreement with reference results from literature. • Robust high-resolution numerical model and algorithm for compressible multi-phase flow with arbitrary number of phases. • No further time-step restrictions implied by high-order interface reconstruction. • Accurate capturing and advection of multi-phase interfaces with positivity preservation. • No numerical artifacts for wave propagation across interfaces. [ABSTRACT FROM AUTHOR]

This paper introduces a novel method for the computation of accurate normal fields from volume fractions on unstructured polyhedral meshes. For surfaces of considerable curvature variation, the method is second-order accurate with respect to symmetric volume differences and first-order accurate with respect to normal field angular deviation. For this purpose, an averaged normal is locally computed, i.e. in each mesh cell, by fitting a plane in a least squares sense with appropriate weights to the volume fraction data of neighbouring cells while explicitly accounting for volume conservation in the cell at hand. The resulting minimization problem is solved approximately by means of a Newton -type method. Moreover, employing the Reynolds transport theorem allows for both assessing the regularity of the derivatives of the least squares error functional with respect to the free parameters and to exactly compute these derivatives at small additional costs. Since the divergence theorem implies that the volume fraction can be written as a sum of face-based quantities, our method considerably simplifies the numerical procedure for applications in three spatial dimensions while demonstrating an inherent ability to robustly deal with unstructured meshes. We discuss the theoretical foundations, regularity and appropriate error measures, along with the details of the numerical algorithm. Finally, numerical results for convex and non-convex hypersurfaces embedded in cuboidal and tetrahedral meshes are presented together with a convergence study considering two different error measures. Moreover, our findings are complemented by profound insights into the minimization procedure. • Efficient second-order accurate computation of normal fields solely from volume fractions. • Applicable to unstructured grids with arbitrary polyhedral cells in three spatial dimensions. • Face-based normal reconstruction with intrinsic exact volume conservation. • Dynamic stencil extension and bulk cell weighting for removal of outliers. • Second-order convergence in symmetric volume difference for moderate curvature variations. [ABSTRACT FROM AUTHOR]

Partial discharges and the electrical treeing represent one of key aspects for the ageing of the insulating parts of electrical components. To predict the end of life of dielectric systems several different modelling approaches have been developed. One of the most promising is the modelling using first principles and a set of partial differential equations. In this field, one of the toughest challenges is to cope with the very different scales of the treeing (each branch is a few micrometers wide) and of electrical components (raging from a few centimetres to a few meters). In this paper we propose an approach, based on the finite element and finite volume methods, which uses two completely independent grids to describe internal defects and electrical components. Some interpolation algorithms have been developed to couple the solution on the two grids. We have also tested the accuracy of our approach by comparing it with previously published methods that use only one grid. • Partial discharges simulation using plasma physics for challenging geometries. • Efficient treatment of tiny internal defects in electric equipments. • Grid decoupling between a gaseous defect and a dielectric bulk. [ABSTRACT FROM AUTHOR]

Aircraft aerodynamic design often involves evaluating flow conditions that span low subsonic to transonic or even supersonic Mach numbers. Compressible flow solvers are a natural choice for such design problems, but these solvers suffer from reduced accuracy and efficiency at low Mach numbers. In addition, simulations with supersonic conditions can be challenging to converge because of large gradients in the flow field. This paper presents three contributions to address these issues in the context of an approximate Newton–Krylov solver for the Reynolds-averaged Navier–Stokes equations. First, we propose a method for scaling the artificial dissipation terms in the Jameson–Schmidt–Turkel scheme to improve its accuracy at low Mach numbers while retaining the simplicity of the original scalar dissipation formulation. Second, we show that characteristic time-stepping combined with an approximate Newton method can accelerate convergence for low Mach number flows by reducing the stiffness in the linear system for each Newton iteration. Third, we introduce a dissipation-based continuation method for flows with shocks that improves robustness and accelerates convergence without sacrificing accuracy. These methods can make compressible flow solvers more accurate and efficient across low and high Mach number regimes. • A simple change to a scalar dissipation scheme improves accuracy at low Mach numbers. • Characteristic time-stepping reduces linear system stiffness for Newton-type solvers. • Scaling artificial dissipation is needed to achieve speedups at low Mach numbers. • Dissipation-based continuation improves robustness for flows with shocks. [ABSTRACT FROM AUTHOR]

In this paper, we propose a cusp-capturing physics-informed neural network (PINN) to solve discontinuous-coefficient elliptic interface problems whose solution is continuous but has discontinuous first derivatives on the interface. To find such a solution using neural network representation, we introduce a cusp-enforced level set function as an additional feature input to the network to retain the inherent solution properties; that is, capturing the solution cusps (where the derivatives are discontinuous) sharply. In addition, the proposed neural network has the advantage of being mesh-free, so it can easily handle problems in irregular domains. We train the network using the physics-informed framework in which the loss function comprises the residual of the differential equation together with certain interface and boundary conditions. We conduct a series of numerical experiments to demonstrate the effectiveness of the cusp-capturing technique and the accuracy of the present network model. Numerical results show that even using a one-hidden-layer (shallow) network with a moderate number of neurons and sufficient training data points, the present network model can achieve prediction accuracy comparable with traditional methods. Besides, if the solution is discontinuous across the interface, we can simply incorporate an additional supervised learning task for solution jump approximation into the present network without much difficulty. • Cusp-capturing physics informed neural network for elliptic interface problems. • Our network can present continuous solutions that inherently have discontinuous first derivatives on interfaces. • A mesh-free approach for solving PDEs with interfaces is presented. • Our network can achieve high prediction accuracy with large contrast discontinuous coefficients. [ABSTRACT FROM AUTHOR]

Kernels obtained from the heat equation arise in several modeling contexts, like some double porosity models, or viscous flows in networks of thin tubes. These kernels are weakly singular at initial time. An accurate approximation must therefore take this singularity into account. In this paper we obtain an asymptotic expansion for small times, which we use to build a numerical scheme for approximating the kernels. Convergence of the scheme and relevance of a correction through the asymptotics are proven both analytically and numerically. Finally, we show that this approximation applies to the model on the graph studied by the authors in Canon et al. [4] (2021). • Small time asymptotic expansion of weakly singular kernels • Convergent numerical approximation of a diffusion equation with incompatible initial and boundary data • Application to simplified models for viscous fluids in networks in relation to blood vessels flows and microfluidics [ABSTRACT FROM AUTHOR]