It is well-known in the Smoothed Particle Hydrodynamics (SPH) community that correction in the gradient and Laplacian operators have the potential to drastically increase the accuracy of the method at the expense of computational stability. This paper proposes a stable implementation of such corrections in all derivative operators to the Arbitrary Lagrangian Eulerian incompressible SPH (ALE-ISPH) method, in addition to a novel Neumann boundary condition (BC) applied directly on the velocity (as opposed to traditional BCs where the constraint is applied on the acceleration). In this way, the pressure is solved for both water and wall particles simultaneously, leading to a pressure field that obeys non-penetration BC and divergence-free at the same time. Furthermore, to stabilize the method, we have developed a novel density-based particle shifting technique (PST), specifically designed to deal with incompressible fluids. In this formulation, the numerical density is given as one of the most critical constraint variables. As a result, the proposed density-based PST can maintain the fluid's overall volume for the whole simulation. In addition, it also provides numerical stability as it prevents particle clustering and leads the fluid domain to an isotropic composition. First, we verified the proposed corrected formulation with the novel Neumann BC for both non-penetration and non-slip conditions with the simulation of hydrostatic pressure and Poisenuille flow, respectively. Then, we tested the proposed density-based PST with the rotating square patch problem with results comparable to previous studies. Lastly, we verified the proposed method for the dam break with an obstacle test, a highly dynamic problem.
This paper presents a computational methodology developed for a high-order approximation of compressible fluid dynamics equations with discontinuities. The methodology is based on a discontinuous Galerkin spectral-element method (DGSEM) built upon a split discretization framework with summation-by-parts (SBP) property, which mimics the integration-by-parts operation in a discrete sense. To extend the split DGSEM framework to discontinuous cases, we implement a shock capturing method based on the entropy viscosity formulation. The developed high-order split-form DGSEM with shock-capturing methodology is subject to a series of evaluation on both one-dimensional and two-dimensional, continuous and discontinuous cases. Convergence of the method is demonstrated both for smooth and shocked cases that have analytical solutions. The 2D Riemann problem tests illustrate an accurate representation of all the relevant flow phenomena, such as shocks, contact discontinuities, and rarefaction waves. All test cases are able to run with a polynomial order of 7 or higher. The values of the tunable parameters related to the entropy viscosity are robust for both 1D and 2D test problems. We also show that higher-order approximations yield smaller errors than lower-order approximations, for the same number of total degrees of freedom.
We propose a physics-driven stretching function for direct numerical simulation (DNS) of compressible turbulent wall-bounded flows, which blends uniform near-wall spacing with uniform resolution in terms of semi-local Kolmogorov units in the outer wall layer. Given target Mach number, Reynolds number and wall temperature, our procedure yields a well-defined prescription for the number of grid points and their distribution which guarantee at the same time numerical accuracy and judicious exploitation of computational resources. DNS of high-speed turbulent boundary layers are used to evaluate the quality of the proposed stretching function, which show that one can achieve identical results as with general-purpose stretching functions, however with substantially higher efficiency. A Python script is provided to facilitate implementation of the proposed grid stretching.
We present the three-dimensional version of the Elliptical Parcel-In-Cell (EPIC) method for the simulation of fluid flows and analogous continuum systems. The method represents a flow using a space-filling set of ellipsoidal parcels, which move, rotate and deform in the flow field. Additionally, parcels may carry any number of attributes, such as vorticity, density, temperature, etc, which generally evolve in time on the moving parcels. An underlying grid is used for efficiency in computing the velocity field from the interpolated vorticity field, and in obtaining parcel attribute tendencies. Mixing is enabled by permitting parcels to split when excessively deformed, and by merging very small parcels with the nearest other parcel. Several tests are provided which illustrate the behaviour of the method and demonstrate its effectiveness in modelling complex, buoyancy-driven turbulent fluid flows. The results are compared with a large eddy simulation (LES) and a direct numerical simulation (DNS) model.
Geometric conduction blocks stop cardiac electric propagation due to the shape or conductivity properties of the domain. The blocks are considered to cause many abnormal cardiac electric propagations, leading to cardiac electrophysiological pathologies, such as cardiac fibrillation and arrhythmia. Locating such multidimensional conduction blocks is challenging, particularly in a complex domain with a complex shape and strong anisotropy, such as the heart. To address this problem, we propose a novel mathematical model of the geometric conduction block using the relative acceleration adopted from space-time physics. An efficient numerical scheme for the mathematical model is also proposed to predict the unidirectional conduction block effectively, even in a complex domain. The relative acceleration in the cardiac electric propagation corresponds to the sink-source relationship between the excited (after repolarization) and excitable (before depolarization) cardiac cells, representing the geometric growth rate of the volume of metric balls. The trajectory is constructed from the wavefront of diffusion-reaction equations by aligning orthonormal basis vectors along the gradient of the action potential. Relative acceleration is computed along the propagational direction from the connection 1-form of the basis vectors. The proposed mathematical model and numerical scheme are applied to demonstrate geometric conduction blocks in two-dimensional (2D) simple curved domains with strong anisotropy.
We present an extension to the well-known peel off optimization for Monte Carlo radiative transfer simulations. The classical method is only applicable when the distance between the detector and the peel off event is much bigger than the size of the detector. We use two alternatives to the classical method and calculate the peel off intensity via a subdivision and an integration method. We compare their performance for a realistic scenario and derive guidelines for a general treatment. This allows for precise peel off calculations at any distance to the detector surface.
Herein we describe a new approach to modelling inviscid two-dimensional stratified flows in a general domain. The approach makes use of a conformal map of the domain to a rectangle. In this transformed domain, the equations of motion are largely unaltered, and in particular Laplace's equation remains unchanged. This enables one to construct exact solutions to Laplace's equation and thereby enforce all boundary conditions.
The nonlinear Heisenberg-Euler theory is capable of describing the dynamics of vacuum polarization, a key prediction by quantum electrodynamics. Due to vast progress in the field of laser technology in recent years vacuum polarization can be triggered in the lab by colliding high-intensity laser pulses, leading to a variety of interesting novel phenomena. Since analytical methods for highly nonlinear problems are generally limited and since the experimental requirements for the detection of the signals from the nonlinear quantum vacuum are high, the need for numerical support is apparent. The paper presents a highly-accurate, efficient numerical scheme for solving the nonlinear Heisenberg-Euler equations in weak-field expansion up to six-photon interactions. Properties of the numerical scheme are discussed and an implementation accurate up to order thirteen in terms of spatial resolution is given. Simulations are presented and benchmarked with known analytical results. The versatility of the numerical solver is demonstrated by solving problems in complicated configurations.
We present a multi-physics model for the approximation of the coupled system formed by the heat equation and the Navier-Stokes equations with solidification and free surfaces. The computational domain is the union of two overlapping regions: a larger domain to account for thermal effects, and a smaller region to account for the fluid flow. Temperature-dependent surface effects are accounted for via surface tension and Marangoni forces. The volume-of-fluid approach is used to track the free surfaces between the metal (liquid or solidified) and the ambient air. The numerical method incorporates all the physical phenomena within an operator splitting strategy. The discretization relies on a two-grid approach that uses an unstructured finite element mesh for diffusion phenomena and a structured Cartesian grid for advection phenomena. The model is validated through numerical experiments, the main application being laser melting and polishing.
Smoothing is a specialized form of Bayesian inference for state-space models that characterizes the posterior distribution of a collection of states given an associated sequence of observations. Ramgraber et al. (2023) proposes a general framework for transport-based ensemble smoothing, which includes linear Kalman-type smoothers as special cases. Here, we build on this foundation to realize and demonstrate nonlinear backward ensemble transport smoothers. We discuss parameterization and regularization of the associated transport maps, and then examine the performance of these smoothers for nonlinear and chaotic dynamical systems that exhibit non-Gaussian behavior. In these settings, our nonlinear transport smoothers yield lower estimation error than conventional linear smoothers and state-of-the-art iterative ensemble Kalman smoothers, for comparable numbers of model evaluations.
Smoothers are algorithms for Bayesian time series re-analysis. Most operational smoothers rely either on affine Kalman-type transformations or on sequential importance sampling. These strategies occupy opposite ends of a spectrum that trades computational efficiency and scalability for statistical generality and consistency: non-Gaussianity renders affine Kalman updates inconsistent with the true Bayesian solution, while the ensemble size required for successful importance sampling can be prohibitive. This paper revisits the smoothing problem from the perspective of measure transport, which offers the prospect of consistent prior-to-posterior transformations for Bayesian inference. We leverage this capacity by proposing a general ensemble framework for transport-based smoothing. Within this framework, we derive a comprehensive set of smoothing recursions based on nonlinear transport maps and detail how they exploit the structure of state-space models in fully non-Gaussian settings. We also describe how many standard Kalman-type smoothing algorithms emerge as special cases of our framework. A companion paper (Ramgraber et al., 2023) explores the implementation of nonlinear ensemble transport smoothers in greater depth.
We present an algorithm for compressing the radiosity view factor model commonly used in radiation heat transfer and computer graphics. We use a format inspired by the hierarchical off-diagonal low rank format, where elements are recursively partitioned using a quadtree or octree and blocks are compressed using a sparse singular value decomposition -- the hierarchical matrix is assembled using dynamic programming. The motivating application is time-dependent thermal modeling on vast planetary surfaces, with a focus on permanently shadowed craters which receive energy through indirect irradiance. In this setting, shape models are comprised of a large number of triangular facets which conform to a rough surface. At each time step, a quadratic number of triangle-to-triangle scattered fluxes must be summed; that is, as the sun moves through the sky, we must solve the same view factor system of equations for a potentially unlimited number of time-varying righthand sides. We first conduct numerical experiments with a synthetic spherical cap-shaped crater, where the equilibrium temperature is analytically available. We also test our implementation with triangle meshes of planetary surfaces derived from digital elevation models recovered by orbiting spacecrafts. Our results indicate that the compressed view factor matrix can be assembled in quadratic time, which is comparable to the time it takes to assemble the full view matrix itself. Memory requirements during assembly are reduced by a large factor. Finally, for a range of compression tolerances, the size of the compressed view factor matrix and the speed of the resulting matrix vector product both scale linearly (as opposed to quadratically for the full matrix), resulting in orders of magnitude savings in processing time and memory space.
A fast and reliable geometry optimization algorithm is presented that optimizes atomic positions and lattice vectors simultaneously. Using a series of benchmarks, it is shown that the method presented in this paper outperforms in most cases the standard optimization methods implemented in popular codes such as Quantum ESPRESSO and VASP. To motivate the variable cell shape optimization method presented in here, the eigenvalues of the lattice Hessian matrix are investigated thoroughly. It is shown that they change depending on the shape of the cell and the number of particles inside the cell. For certain cell shapes the resulting condition number of the lattice matrix can grow quadratically with respect to the number of particles. By a coordinate transformation, which can be applied to all variable cell shape optimization methods, the undesirable conditioning of the lattice Hessian matrix is eliminated.
Mechanochemical processes on surfaces such as the cellular cortex or epithelial sheets, play a key role in determining patterns and shape changes of biological systems. To understand the complex interplay of hydrodynamics and material flows on such active surfaces requires novel numerical tools. Here, we present a finite-element method for an active deformable surface interacting with the surrounding fluids. The underlying model couples surface and bulk hydrodynamics to surface flow of a diffusible species which generates active contractile forces. The method is validated with previous results based on linear stability analysis and shows almost perfect agreement regarding predicted patterning. Away from the linear regime we find rich non-linear behavior, such as the presence of multiple stationary states. We study the formation of a contractile ring on the surface and the corresponding shape changes. Finally, we explore mechanochemical pattern formation on various surface geometries and find that patterning strongly adapts to local surface curvature. The developed method provides a basis to analyze a variety of systems that involve mechanochemical pattern formation on active surfaces interacting with surrounding fluids.
Numerical models for predicting future ice mass loss of the Antarctic and Greenland ice sheets require accurately representing their dynamics. Unfortunately, ice-sheet models suffer from a very strict time-step size constraint, which for higher-order models constitutes a severe bottleneck; in each time step a nonlinear and computationally demanding system of equations has to be solved. In this study, stable time-step sizes are increased for a full-Stokes model by implementing a so-called free-surface stabilization algorithm (FSSA). Previously this stabilization has been used successfully in mantle-convection simulations where a similar viscous-flow problem is solved. By numerical investigation it is demonstrated that instabilities on the very thin domains required for ice-sheet modeling behave differently than on the equal-aspect-ratio domains the stabilization has previously been used on. Despite this, and despite the different material properties of ice, it is shown that it is possible to adapt FSSA to work on idealized ice-sheet domains and increase stable time-step sizes by at least one order of magnitude. The FSSA method presented is deemed accurate, efficient and straightforward to implement into existing ice-sheet solvers.
Estimating modeling parameters based on a prescribed optimization target requires to solve an inverse problem, which is commonly ill-posed. Consequently, either infinitely many or no solutions may exist, depending on whether the system is under- or overdetermined, and whether it is consistent or inconsistent. This paper focuses on scenarios where the solution is ambiguous and infinitely many combinations of possible parameter values can accurately achieve the optimization target. Selecting the most suitable solution requires incorporating additional constraints into the model, which is achieved by regularizing the inverse problem. However, common regularization approaches require the specification of a priori unknown regularization hyperparameters that are difficult and tedious to obtain, and can have a large impact on the result. Here, a novel strategy to reduce the ambiguity of such inverse problems is presented, ensuring that the primary optimization target is always reached accurately. To further reduce the solution space, additional constraints are included, until the optimal modeling parameters are found. Importantly, the required regularization parameters have a direct physical meaning and can be derived sequentially, starting from an initial guess that can be obtained conveniently by solving the system without regularization. By considering several illustrative examples, the applicability of the method is demonstrated, and its potential for various comparable inverse problems is highlighted.
We present a novel approach to simulating general two-dimensional flows, which could also be applied to other areas of continuum mechanics. The approach generalises the Particle-In-Cell (PIC) method, originally used to model two-dimensional hydrodynamics, by representing fluid elements by elliptical parcels. The rotation and deformation of these parcels are calculated, and parcels split beyond a critical aspect ratio. Conversely, small parcels are eliminated by merging them with larger ones. The elliptical parcels well represent the flow deformation and have excellent conservation properties. In contrast to earlier work that combined PIC with elliptical parcels that split and merge, a vorticity-based framework is used, and accurate integration over ellipses is performed efficiently by two-point Gaussian quadrature. The small-scale mixing associated with parcel splitting and merging is shown to be strongly convergent with grid resolution. The robustness, versatility, accuracy and efficiency of the new Elliptical Parcel-In-Cell (EPIC) method is demonstrated for a variety of standard test cases, and compared with a standard pseudo-spectral method. The results indicate that EPIC is a promising, Lagrangian-based alternative to grid-based methods.
This paper presents a spectral scheme for the numerical solution of nonlinear conservation laws in non-periodic domains under arbitrary boundary conditions. The approach relies on the use of the Fourier Continuation (FC) method for spectral representation of non-periodic functions in conjunction with smooth localized artificial viscosity assignments produced by means of a Shock-Detecting Neural Network (SDNN). Like previous shock capturing schemes and artificial viscosity techniques, the combined FC-SDNN strategy effectively controls spurious oscillations in the proximity of discontinuities. Thanks to its use of a localized but smooth artificial viscosity term, whose support is restricted to a vicinity of flow-discontinuity points, the algorithm enjoys spectral accuracy and low dissipation away from flow discontinuities, and, in such regions, it produces smooth numerical solutions -- as evidenced by an essential absence of spurious oscillations in level set lines. The FC-SDNN viscosity assignment, which does not require use of problem-dependent algorithmic parameters, induces a significantly lower overall dissipation than other methods, including the Fourier-spectral versions of the previous entropy viscosity method. The character of the proposed algorithm is illustrated with a variety of numerical results for the linear advection, Burgers and Euler equations in one and two-dimensional non-periodic spatial domains.
Scientific machine learning has been successfully applied to inverse problems and PDE discovery in computational physics. One caveat concerning current methods is the need for large amounts of ("clean") data, in order to characterize the full system response and discover underlying physical models. Bayesian methods may be particularly promising for overcoming these challenges, as they are naturally less sensitive to the negative effects of sparse and noisy data. In this paper, we propose to use Bayesian neural networks (BNN) in order to: 1) Recover the full system states from measurement data (e.g. temperature, velocity field, etc.). We use Hamiltonian Monte-Carlo to sample the posterior distribution of a deep and dense BNN, and show that it is possible to accurately capture physics of varying complexity, without overfitting. 2) Recover the parameters instantiating the underlying partial differential equation (PDE) governing the physical system. Using the trained BNN, as a surrogate of the system response, we generate datasets of derivatives that are potentially comprising the latent PDE governing the observed system and then perform a sequential threshold Bayesian linear regression (STBLR), between the successive derivatives in space and time, to recover the original PDE parameters. We take advantage of the confidence intervals within the BNN outputs, and introduce the spatial derivatives cumulative variance into the STBLR likelihood, to mitigate the influence of highly uncertain derivative data points; thus allowing for more accurate parameter discovery. We demonstrate our approach on a handful of example, in applied physics and non-linear dynamics.
While fast multipole methods (FMMs) are in widespread use for the rapid evaluation of potential fields governed by the Laplace, Helmholtz, Maxwell or Stokes equations, their coupling to high-order quadratures for evaluating layer potentials is still an area of active research. In three dimensions, a number of issues need to be addressed, including the specification of the surface as the union of high-order patches, the incorporation of accurate quadrature rules for integrating singular or weakly singular Green's functions on such patches, and their coupling to the oct-tree data structures on which the FMM separates near and far field interactions. Although the latter is straightforward for point distributions, the near field for a patch is determined by its physical dimensions, not the distribution of discretization points on the surface. Here, we present a general framework for efficiently coupling locally corrected quadratures with FMMs, relying primarily on what are called generalized Gaussian quadratures rules, supplemented by adaptive integration. The approach, however, is quite general and easily applicable to other schemes, such as Quadrature by Expansion (QBX). We also introduce a number of accelerations to reduce the cost of quadrature generation itself, and present several numerical examples of acoustic scattering that demonstrate the accuracy, robustness, and computational efficiency of the scheme. On a single core of an Intel i5 2.3GHz processor, a Fortran implementation of the scheme can generate near field quadrature corrections for between 1000 and 10,000 points per second, depending on the order of accuracy and the desired precision. A Fortran implementation of the algorithm described in this work is available at https://gitlab.com/fastalgorithms/fmm3dbie.
