Journal of Computational Physics: X最新文献

Source localization based on electroencephalography (EEG) has become a widely used neuroimaging technique. However its precision has been shown to be very dependent on how accurately the brain, head and scalp can be electrically modeled within the so-called forward problem. The construction of this model is traditionally performed by leveraging Finite Element or Boundary Element Methods (FEM or BEM). Even though the latter is more computationally efficient thanks to the smaller interaction matrices it yields and near-linear solvers, it has traditionally been used on simpler models than the former. Indeed, while FEM models taking into account the different media anisotropies are widely available, BEM models have been limited to isotropic, piecewise homogeneous models. In this work we augment standard BEM with a new wire integral equation to account for the anisotropy of the white matter. The new formulation combines the efficiency of BEM discretization of the boundaries only and modeling of the fibrous nature of the white matter as one-dimensional basis functions which limits the computational impact of their modeling. We compare our scheme against widely used formulations and establish its correctness in both canonical and realistic cases.

We present an adaptive local time-stepping (ALTS) scheme for a block-structured multiresolution scheme of hyperbolic conservation laws for fluid flow. The stability of standard local time-stepping (LTS) schemes with level-dependent time-step sizes is improved by local time-step size adaptation when progressing through the underlying multi-stage time integration scheme. The novelty of the approach is that it merges flux computation and time integration of the state vector with projection and prediction operations of the multiresolution scheme [15]. This enables consistent time integration of subdomains with different refinement levels without the need for intermediate time synchronization which can be prohibitively expensive in parallel computations. Consequently, coarser subdomains are advanced in time only once finer subdomains have advanced to the same time instant. Full spatial resolution adaptivity for integrated regions after each substep is maintained.

The new scheme exhibits significantly improved numerical stability as compared to previous LTS schemes due to the local time-step size adaptation at each substep. The computational overhead of the incurred additional operations is small. In applications, the ALTS scheme demonstrates the same computational efficiency as standard LTS schemes.

The new scheme can be applied to any explicit single-step time-integration scheme and is independent of the employed spatial discretization scheme. The improved stability is demonstrated with several one- and two-dimensional examples of flows with one and two phases, applying second- and third-order Runge-Kutta time integration schemes.

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.

The Lorentz equations describe the motion of electrically charged particles in electric and magnetic fields and are used widely in plasma physics. The most popular numerical algorithm for solving them is the Boris method, a variant of the Störmer-Verlet algorithm. Boris method is phase space volume conserving and simulated particles typically remain near the correct trajectory. However, it is only second order accurate. Therefore, in scenarios where it is not enough to know that a particle stays on the right trajectory but one needs to know where on the trajectory the particle is at a given time, Boris method requires very small time steps to deliver accurate phase information, making it computationally expensive. We derive an improved version of the high-order Boris spectral deferred correction algorithm (Boris-SDC) by adopting a convergence acceleration strategy for second order problems based on the Generalised Minimum Residual (GMRES) method. Our new algorithm is easy to implement as it still relies on the standard Boris method. Like Boris-SDC it can deliver arbitrary order of accuracy through simple changes of runtime parameter but possesses better long-term energy stability. We demonstrate for two examples, a magnetic mirror trap and the Solev'ev equilibrium, that the new method can deliver better accuracy at lower computational cost compared to the standard Boris method. While our examples are motivated by tracking ions in the magnetic field of a nuclear fusion reactor, the introduced algorithm can potentially deliver similar improvements in efficiency for other applications.

We construct high-order derivative operators for smooth functions represented via discontinuous multiwavelet bases. The need for such operators arises in order to avoid artifacts when computing functionals involving high-order derivatives of solutions of integral equations. Previously high-order derivatives had to be formed by repeated application of a first-derivative operator that, while uniquely defined, has a spectral norm that grows quadratically with polynomial order and, hence, greatly amplifies numerical noise (truncation error) in the multiwavelet computation. The new constructions proceed via least-squares projection onto smooth bases and provide substantially improved numerical properties as well as permitting direct construction of high-order derivatives. We employ either b-splines or bandlimited exponentials as the intermediate smooth basis, with the former maintaining the concept of approximation order while the latter preserves the pure imaginary spectrum of the first-derivative operator and provides more direct control over the bandlimit and accuracy of computation. We demonstrate the properties of these new operators via several numerical tests as well as application to a problem in nuclear physics.

Active flux schemes are finite volume schemes that keep track of both point values and averages. The point values are updated using a semi-Lagrangian step, making active flux schemes highly suitable for geometric optics problems on phase space, i.e., to solve Liouville's equation. We use a semi-discrete version of the active flux scheme. Curved optics lead to moving boundaries in phase space. Therefore, we introduce a novel way of defining the active flux scheme on moving meshes. We show, using scaling arguments as well as numerical experiments, that our scheme outperforms the current industry standard, ray tracing. It has higher accuracy as well as a more favourable time scaling. The numerical experiments demonstrate that the active flux scheme is orders of magnitude more accurate and faster than ray tracing.

Developing a turbulence model that is computationally inexpensive and compatible with the nature of the numerical scheme is a crucial step in expanding the application of spectral element methods for large eddy simulation (LES) in complex geometries. In this paper, an element-level modal low-pass explicit filtering procedure, which operates in the spectral space, is implemented in a discontinuous spectral element method (DSEM). The application of the modal filter is studied for LES without a subgrid-scale (SGS) model. The method is tested for a configuration featuring isotropic turbulence, and its performance is compared with a previously used method—a spectral interpolation-based nodal filter. The modal filter shows superior performance over the nodal filter. The filtering procedure is then applied to a turbulent channel flow at a friction Reynolds number of ${\mathrm{Re}}_{\tau }=544$, and the results are compared with a previous direct numerical simulation (DNS). It is also shown that the filter strength that provides the best comparison with DNS depends only on the polynomial order and is not a function of the grid resolution. An anisotropic version of the modal filter, which damps high-frequency modes in a specific direction, is also introduced and tested for the channel flow. It is observed that filtering in the spanwise direction is the most effective approach based on the comparison of velocity mean and fluctuations with DNS. In general, the modal filter has shown good performance for both isotropic and wall-bounded flows; the calculated channel friction Reynolds number for the modal filter is within 0.26% error with respect to the DNS data, compared to 5.8% error for a case with no filtering.

This paper presents a method for accurate and efficient computations on scalar, vector and tensor fields in three-dimensional spherical polar coordinates. The method uses spin-weighted spherical harmonics in the angular directions and rescaled Jacobi polynomials in the radial direction. For the 2-sphere, spin-weighted harmonics allow for automating calculations in a fashion as similar to Fourier series as possible. Derivative operators act as wavenumber multiplication on a set of spectral coefficients. After transforming the angular directions, a set of orthogonal tensor rotations put the radially dependent spectral coefficients into individual spaces each obeying a particular regularity condition at the origin. These regularity spaces have remarkably simple properties under standard vector-calculus operations, such as gradient and divergence. We use a hierarchy of rescaled Jacobi polynomials for a basis on these regularity spaces. It is possible to select the Jacobi-polynomial parameters such that all relevant operators act in a minimally banded way. Altogether, the geometric structure allows for the accurate and efficient solution of general partial differential equations in the unit ball.

This paper presents an explicit vertex centred finite volume method for the solution of fast transient isothermal large strain solid dynamics via a system of first order hyperbolic conservation laws. Building upon previous work developed by the authors, in the context of alternative numerical discretisations, this paper explores the use of a series of enhancements (both from the formulation and numerical standpoints) in order to explore some limiting scenarios, such as the consideration of near and true incompressibility. Both Total and Updated Lagrangian formulations are presented and compared at the discrete level, where very small differences between both descriptions are observed due to the excellent discrete satisfaction of the involutions. In addition, a matrix-free tailor-made artificial compressibility algorithm is discussed and combined with an angular momentum projection algorithm. A wide spectrum of numerical examples is thoroughly examined. The scheme shows excellent (stable, consistent and accurate) behaviour, in comparison with other methodologies, in compressible, nearly incompressible and truly incompressible bending dominated scenarios, yielding equal second order of convergence for velocities, deviatoric and volumetric components of the stress.

We present a simulation code which can solve a broad range of partial differential equations in a full sphere. The code expands tensorial variables in a spectral series of spin-weighted spherical harmonics in the angular directions and a scaled Jacobi polynomial basis in the radial direction, as described in Vasil et al. (2018; hereafter, Part-I). Nonlinear terms are calculated by transforming from the coefficients of the spectral series to the value of each quantity on the physical grid, where it is easy to calculate products and perform other local operations. The expansion makes it straightforward to solve equations in tensor form (i.e., without decomposition into scalars). We propose and study several unit tests which demonstrate the code can accurately solve linear problems, implement boundary conditions, and transform between spectral and physical space. We then run a series of benchmark problems proposed in Marti et al. (2014), implementing the hydrodynamic and magnetohydrodynamic equations. We are able to calculate more accurate solutions than reported in Marti et al. (2014) by running at higher spatial resolution and using a higher-order timestepping scheme. We find the rotating convection and convective dynamo benchmark problems depend sensitively on details of timestepping and data analysis. We also demonstrate that in low resolution simulations of the dynamo problem, small changes in a numerical scheme can lead to large changes in the solution. To aid future comparison to these benchmarks, we include the source code used to generate the data, as well as the data and analysis scripts used to generate the figures.