ALBERT

Description: 〈p〉Publication date: Available online 8 July 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Maxim Rakhuba, Alexander Novikov, Ivan Oseledets〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Such problems as computation of spectra of spin chains and vibrational spectra of molecules can be written as 〈em〉high-dimensional eigenvalue problems〈/em〉, i.e., when the eigenvector can be naturally represented as a multidimensional tensor. Tensor methods have proven to be an efficient tool for the approximation of solutions of high-dimensional eigenvalue problems, however, their performance deteriorates quickly when the number of eigenstates to be computed increases. We address this issue by designing a new algorithm motivated by the ideas of 〈em〉Riemannian optimization〈/em〉 (optimization on smooth manifolds) for the approximation of multiple eigenstates in the 〈em〉tensor-train format〈/em〉, which is also known as matrix product state representation. The proposed algorithm is implemented in TensorFlow, which allows for both CPU and GPU parallelization.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 8 July 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Chen Liu, Florian Frank, Faruk O. Alpak, Béatrice Rivière〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Permeability estimation of porous media from directly solving the Navier–Stokes equations has a wide spectrum of applications in petroleum industry. In this paper, we utilize a pressure-correction projection algorithm in conjunction with the interior penalty discontinuous Galerkin scheme for space discretization to build an incompressible Navier–Stokes simulator and to use this simulator to calculate permeability of real rock samples. The proposed method is accurate, numerically robust, and exhibits the potential for tackling realistic problems.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 8 July 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Mustapha Malek, Nouh Izem, M. Shadi Mohamed, Mohammed Seaid, Omar Laghrouche〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉An efficient partition of unity finite element method for three-dimensional transient diffusion problems is presented. A class of multiple exponential functions independent of time variable is proposed to enrich the finite element approximations. As a consequence of this procedure, the associated matrix for the linear system is evaluated once at the first time step and the solution is obtained at subsequent time step by only updating the right-hand side of the linear system. This results in an efficient numerical solver for transient diffusion equations in three space dimensions. Compared to the conventional finite element methods with 〈em〉h〈/em〉-refinement, the proposed approach is simple, more efficient and more accurate. The performance of the proposed method is assessed using several test examples for transient diffusion in three space dimensions. We present numerical results for a transient diffusion equation with known analytical solution to quantify errors for the new method. We also solve time-dependent diffusion problems in complex geometries. We compare the results obtained using the partition of unity finite element method to those obtained using the standard finite element method. It is shown that the proposed method strongly reduces the necessary number of degrees of freedom to achieve a prescribed accuracy.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 8 July 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Lahbib Bourhrara〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This document presents a new numerical scheme dealing with the Boltzmann transport equation. This scheme is based on the expansion of the angular flux in a truncated spherical harmonics function and the discontinuous finite element method for the spatial variable. The advantage of this scheme lies in the fact that we can deal with unstructured, non-conformal and curved meshes. Indeed, it is possible to deal with distorted regions whose boundary is constituted by edges that can be either line segments or circular arcs or circles. In this document, we detail the derivation of the method for 2D geometries. However, the generalization to 2D extruded geometries is trivial.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 396〈/p〉 〈p〉Author(s): Luigi Brugnano, Juan I. Montijano, Luis Rández〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this paper we study arbitrarily high-order energy-conserving methods for simulating the dynamics of a charged particle. They are derived and studied within the framework of 〈em〉Line Integral Methods (LIMs)〈/em〉, previously used for defining 〈em〉Hamiltonian Boundary Value Methods (HBVMs)〈/em〉, a class of energy-conserving Runge-Kutta methods for Hamiltonian problems. A complete analysis of the new methods is provided, which is confirmed by a few numerical tests.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 21 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Konstantin Pieper, K. Chad Sockwell, Max Gunzburger〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A framework for exponential time discretization of the multilayer rotating shallow water equations is developed in combination with a mimetic discretization in space. The method is based on a combination of existing exponential time differencing (ETD) methods and a careful choice of approximate Jacobians. The discrete Hamiltonian structure and conservation properties of the model are taken into account, in order to ensure stability of the method for large time steps and simulation horizons. In the case of many layers, further efficiency can be gained by a layer reduction which is based on the vertical structure of fast and slow modes. Numerical experiments on the example of a mid-latitude regional ocean model confirm long term stability for time steps increased by an order of magnitude over the explicit CFL, while maintaining accuracy for key statistical quantities.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 20 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Stéphane Zaleski, Feng Xiao〈/p〉

Description: 〈p〉Publication date: Available online 20 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Dario Collia, Marija Vukicevic, Valentina Meschini, Luigino Zovatto, Gianni Pedrizzetti〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉The fluid dynamics inside the left ventricle of the human heart is considered a potential indicator of long term cardiovascular outcome. In this respect, numerical simulations can play an important role for integrating existing technology to reproduce flow details and even conditions associated to virtual therapeutic solutions. Nevertheless, numerical models encounter serious practical difficulties in describing the interaction between flow and surrounding tissues due to the limited information inherently available in real clinical applications.〈/p〉 〈p〉This study presents a computational method for the fluid dynamics inside the left ventricle designed to be efficiently integrated in clinical scenarios. It includes an original model of the mitral valve dynamics, which describes an asymptotic behavior for tissues with no elastic stiffness other than the constrain of the geometry obtained from medical imaging; in particular, the model provides an asymptotic description without requiring details of tissue properties that may not be measurable 〈em〉in vivo〈/em〉.〈/p〉 〈p〉The advantages of this model with respect to a valveless orifice and its limitations with respect to a complete tissue modeling are verified. Its performances are then analyzed in details to ensure a correct interpretation of results. It represents a potential option when information about tissue mechanical properties is insufficient for the implementations of a full fluid-structure interaction approach.〈/p〉 〈/div〉

Description: 〈p〉Publication date: Available online 22 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): L. Nouveau, M. Ricchiuto, G. Scovazzi〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We propose an extension of the embedded boundary method known as “shifted boundary method” to elliptic diffusion equations in mixed form (e.g., Darcy flow, heat diffusion problems with rough coefficients, etc.). Our aim is to obtain an improved formulation that, for linear finite elements, is at least second-order accurate for both flux and primary variable, when either Dirichlet or Neumann boundary conditions are applied. Following previous work of Nishikawa and Mazaheri in the context of residual distribution methods, we consider the mixed form of the diffusion equation (i.e., with Darcy-type operators), and introduce an enrichment of the primary variable. This enrichment is obtained exploiting the relation between the primary variable and the flux variable, which is explicitly available at nodes in the mixed formulation. The proposed enrichment mimics a formally quadratic pressure approximation, although only nodal unknowns are stored, similar to a linear finite element approximation. We consider both continuous and discontinuous finite element approximations and present two approaches: a non-symmetric enrichment, which, as in the original references, only improves the consistency of the overall method; and a symmetric enrichment, which enables a full error analysis in the classical finite element context. Combined with the shifted boundary method, these two approaches are extended to high-order embedded computations, and enable the approximation of both primary and flux (gradient) variables with second-order accuracy, independently on the type of boundary conditions applied. We also show that the the primary variable is third-order accurate, when pure Dirichlet boundary conditions are embedded.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 21 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Walter Boscheri, Dinshaw S. Balsara〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉In this work we present a conservative WENO Adaptive Order (AO) reconstruction operator applied to an explicit one-step Arbitrary-Lagrangian-Eulerian (ALE) discontinuous Galerkin (DG) method. The spatial order of accuracy is improved by reconstructing higher order piecewise polynomials of degree 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈mi〉M〈/mi〉〈mo linebreak="goodbreak" linebreakstyle="after"〉〉〈/mo〉〈mi〉N〈/mi〉〈/math〉, starting from the underlying polynomial solution of degree 〈em〉N〈/em〉 provided by the DG scheme. High order of accuracy in time is achieved by the ADER approach, making use of an element-local space-time Galerkin finite element predictor that arises from a one-step time integration procedure. As a result, space-time polynomials of order 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.svg"〉〈mi〉M〈/mi〉〈mo linebreak="goodbreak" linebreakstyle="after"〉+〈/mo〉〈mn〉1〈/mn〉〈/math〉 are obtained and used to perform the time evolution of the numerical solution adopting a fully explicit DG scheme.〈/p〉 〈p〉To maintain algorithm simplicity, the mesh motion is restricted to be carried out using straight lines, hence the old mesh configuration at time 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si3.svg"〉〈msup〉〈mrow〉〈mi〉t〈/mi〉〈/mrow〉〈mrow〉〈mi〉n〈/mi〉〈/mrow〉〈/msup〉〈/math〉 is connected with the new one at time 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si4.svg"〉〈msup〉〈mrow〉〈mi〉t〈/mi〉〈/mrow〉〈mrow〉〈mi〉n〈/mi〉〈mo linebreak="badbreak" linebreakstyle="after"〉+〈/mo〉〈mn〉1〈/mn〉〈/mrow〉〈/msup〉〈/math〉 via space-time segments, which result in space-time control volumes on which the governing equations have to be integrated in order to obtain the time evolution of the discrete solution. Our algorithm falls into the category of 〈em〉direct〈/em〉 Arbitrary-Lagrangian-Eulerian (ALE) schemes, where the governing PDE system is directly discretized relying on a space-time conservation formulation and which already takes into account the new grid geometry 〈em〉directly〈/em〉 during the computation of the numerical fluxes. A local rezoning strategy might be used in order to locally optimize the mesh quality and avoiding the generation of invalid elements with negative determinant. The proposed approach reduces to direct ALE finite volume schemes if 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si5.svg"〉〈mi〉N〈/mi〉〈mo linebreak="goodbreak" linebreakstyle="after"〉=〈/mo〉〈mn〉0〈/mn〉〈/math〉, while explicit direct ALE DG schemes are recovered in the case of 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si6.svg"〉〈mi〉N〈/mi〉〈mo linebreak="goodbreak" linebreakstyle="after"〉=〈/mo〉〈mi〉M〈/mi〉〈/math〉.〈/p〉 〈p〉In order to stabilize the DG solution, an 〈em〉a priori〈/em〉 WENO based limiting technique is employed, that makes use of the numerical solution inside the element under consideration and its neighbor cells to find a less oscillatory polynomial approximation. By using a 〈em〉modal basis〈/em〉 in a reference element, the evaluation of the oscillation indicators is very easily and efficiently carried out, hence allowing higher order modes to be properly limited, while leaving the zero-〈em〉th〈/em〉 order mode untouched for ensuring conservation.〈/p〉 〈p〉Numerical convergence rates for 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si7.svg"〉〈mn〉2〈/mn〉〈mo〉≤〈/mo〉〈mi〉N〈/mi〉〈mo〉,〈/mo〉〈mi〉M〈/mi〉〈mo〉≤〈/mo〉〈mn〉4〈/mn〉〈/math〉 are presented as well as a wide set of benchmark test problems for hydrodynamics on moving and fixed unstructured meshes.〈/p〉 〈/div〉

Description: 〈p〉Publication date: Available online 20 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Ali Zidane, Abbas Firoozabadi〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Non-planar fractures are often created in hydraulic fracturing. These irregular shape fractures may reduce the penetration into the formation; they may also improve the reservoir reach. Accurate flow simulation of two-phase compositional flows in domains with complex non-planar fractures is beyond the capabilities of current numerical models. In this work we present a higher-order numerical model for compositional two-phase flow in a domain with non-planar fractures. Fully unstructured gridding in 3D is a natural choice for description of geometry with irregular fracture shapes. We apply the concept of fracture cross-flow equilibrium (FCFE) in simulations of porous media flows with non-planar fractures. FCFE allows accurate flow and composition calculations at low CPU cost. Our implementation is in the context of the hybridized form of the mass conservative mixed finite element (MFE) and the higher-order discontinuous Galerkin (DG) method. In this work we introduce a simple and effective approach for design of non-planar fractures through the mesh interface that connects computer-aided-design (CAD) software to the mesh generator. In our algorithm we can simulate all ranges of fracture permeability accurately as opposed to other approaches where low permeability fractures affect the accuracy.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Martin Pollack, Michele Pütz, Daniele L. Marchisio, Michael Oevermann, Christian Hasse〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉The evolution of polydisperse systems is governed by population balance equations. A group of efficient solution approaches are the moment methods, which do not solve for the number density function (NDF) directly but rather for a set of its moments. While this is computationally efficient, a specific challenge arises when describing the fluxes across a boundary in phase space for the disappearance of elements, the so-called zero-flux. The main difficulty is the missing NDF-information at the boundary, which most moment methods cannot provide. Relevant physical examples are evaporating droplets, soot oxidation or particle dissolution.〈/p〉 〈p〉In general, this issue can be solved by reconstructing the NDF close to the boundary. However, this was previously only achieved with univariate approaches, i.e. considering only a single internal variable. Many physical problems are insufficiently described by univariate population balance equations, e.g. droplets in sprays often require the temperature or the velocity to be internal coordinates in addition to the size.〈/p〉 〈p〉In this paper, we propose an algorithm, which provides an efficient multivariate approach to calculate the zero-fluxes. The algorithm employs the Extended Quadrature Method of Moments (EQMOM) with Beta and Gamma kernel density functions for the marginal NDF reconstruction and a polynomial or spline for the other conditional dimensions. This combination allows to reconstruct the entire multivariate NDF and based on this, expressions for the disappearance flux are derived. An algorithm is proposed for the whole moment inversion and reconstruction process. It is validated against a suite of test cases with increasing complexity. The influence of the number of kernel density functions and the configuration of the polynomials and splines on the accuracy is discussed. Finally, the associated computational costs are evaluated.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Antoine Vermeil de Conchard, Huina Mao, Romain Rumpler〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Effective treatment of unbounded domains using artificial truncating boundaries are essential in numerical simulation, e.g. using the Finite Element Method (FEM). Among these, Perfectly Matched Layers (PML) have proved to be particularly efficient and flexible. However, an efficient handling of frequency sweeps is not trivial with such absorbing layers since the formulation inherently contains coupled space- and frequency-dependent terms. Using the FEM, this may imply generating system matrices at each step of the frequency sweep. In this paper, an approximation is proposed in order to allow for efficient frequency sweeps. The performance and robustness of the proposed approximation is presented on 2D and 3D acoustic cases. A generic, robust way to truncate the acoustic domain efficiently is also proposed, tested on a range of test cases and for different frequency regions.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Maxime Theillard, David Saintillan〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present a new framework for the efficient simulation of the dynamics of active fluids in complex two- and three-dimensional microfluidic geometries. Focusing on the case of a suspension of microswimmers such as motile bacteria, we adopt a continuum mean-field model based on partial differential equations for the evolution of the concentration, polarization and nematic tensor fields, which are nonlinearly coupled to the Navier-Stokes equations for the fluid flow driven by internal active stresses. A level set method combined with an adaptive mesh refinement scheme on Quad-/Octree grids is used to capture complex domain shapes while refining the solution near boundaries or in the neighborhood of sharp gradients. A hybrid finite volumes/finite differences method is implemented in which the concentration field is treated using finite volumes to ensure mass conservation, while the polarization and nematic alignment fields are treated using a combination of finite differences and finite volumes for enhanced accuracy. The governing equations for these fields are solved along with the Navier-Stokes equations, which are evolved using an unconditionally stable projection solver. We illustrate the versatility and robustness of our method by analyzing spontaneous active flows in various two- and three-dimensional systems. Our results show excellent agreement with previous models and experiments and pave the way for further developments in active microfluidics.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Jiangming Xie, M. Yvonne Ou, Liwei Xu〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉Poroelastic materials play an important role in biomechanical and geophysical research. In this paper, we investigate wave propagation in orthotropic poroelastic media by studying the time-domain poroelastic wave equations. Both the low frequency Biot's (LF-Biot) equations and the Biot-Johnson-Koplik-Dashen (Biot-JKD) model are considered. In LF-Biot equations, the dissipation terms are proportional to the relative velocity between the fluid and the solid by a constant. Contrast to this, the dissipation terms in the Biot-JKD model are in the form of time convolution (memory) as a result of the frequency-dependence of fluid-solid interaction at the underlying microscopic scale in the frequency domain. The dynamic tortuosity and permeability described by Darcy's law are two crucial factors in this problem, and highly linked to the viscous force. In the Biot model, the key difficulty is to handle the viscous term when the pore fluid is viscous. In the Biot-JKD model, the convolution operator involves order 1/2 shifted fractional derivatives in the time domain, which is challenging to discretize.〈/p〉 〈p〉In this work, a new method of the multipoint Padé (or Rational) approximation for Stieltjes function is applied to approximate the JKD dynamic tortuosity and then an augmented system of Biot-JKD model is obtained, where the kernel of the memory term is replaced by the finite auxiliary variables satisfying a local system of ordinary differential equations. The Runge-Kutta discontinuous Galerkin (RKDG) method with the un-splitting method is used to compute the numerical solution, and numerical examples are presented to demonstrate the high order accuracy and stability of the method. Compared with the existing approaches for solving the Biot-JKD equations, the augmented system presented here require neither the storage of solution history nor the computation of the flux of the auxiliary variables.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): William C. Tyson, Gary K. Yan, Christopher J. Roy, Carl F. Ollivier-Gooch〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A higher-order accurate discretization error estimation procedure for finite-volume schemes is presented. Discretization error estimates are computed using the linearized error transport equations (ETE). ETE error estimates are applied as a correction to the primal solution. The ETE are then relinearized about the corrected primal solution, and discretization error estimates are recomputed. This process, referred to as ETE relinearization, is performed in an iterative manner to successively increase the accuracy of discretization error estimates. Under certain conditions, ETE relinearization is shown to correct error estimates, or equivalently the entire primal solution, to higher-order accuracy. In terms of computational cost, ETE relinearization has a competitive advantage over conventional higher-order discretizations when used as a form of defect correction for the primal solution. Furthermore, ETE relinearization is shown to be particularly useful for problems where the error incurred by the linearization of the ETE cannot be neglected. Results are presented for several steady-state inviscid and viscous flow problems using both structured and unstructured meshes.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Yu Li, Richard Mikaël Slevinsky〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present a unified treatment of the Fourier spectra of spherically symmetric nonlocal diffusion operators. We develop numerical and analytical results for the class of kernels with weak algebraic singularity as the distance between source and target tends to 0. Rapid algorithms are derived for their Fourier spectra with the computation of each eigenvalue independent of all others. The algorithms are trivially parallelizable, capable of leveraging more powerful compute environments, and the accuracy of the eigenvalues is individually controllable. The algorithms include a Maclaurin series and a full divergent asymptotic series valid for any 〈em〉d〈/em〉 spatial dimensions. Using Drummond's sequence transformation, we prove linear complexity recurrence relations for degree-graded sequences of numerators and denominators in the rational approximations to the divergent asymptotic series. These relations are important to ensure that the algorithms are efficient, and also increase the numerical stability compared with the conventional algorithm with quadratic complexity.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Hua Shen, Matteo Parsani〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We construct a space-time conservation element and solution element (CESE) scheme for solving the compressible Euler equations on moving meshes (CESE-MM) which allow an arbitrary motion for each of the mesh points. The scheme is a direct extension of a purely Eulerian CESE scheme that was previously implemented on hybrid unstructured meshes (Shen et al. (2015) [43]). It adopts a staggered mesh in space and time such that the physical variables are continuous across the interfaces of the adjacent space-time control volumes and, therefore, a Riemann solver is not required to calculate interface fluxes or the node velocities. Moreover, the staggered mesh can significantly alleviate mesh tangles so that the time step can be kept at an acceptable level without using any rezoning operation. The discretization of the integral space-time conservation law is completely based on the physical space-time control volume, thereby satisfying the physical and geometrical conservation laws. Plenty of numerical examples are carried out to validate the accuracy and robustness of the CESE-MM scheme.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): S. Dargaville, A.G. Buchan, R.P. Smedley-Stevenson, P.N. Smith, C.C. Pain〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This paper describes an angular adaptivity algorithm for Boltzmann transport applications which uses P〈sub〉〈em〉n〈/em〉〈/sub〉 and filtered P〈sub〉〈em〉n〈/em〉〈/sub〉 expansions, allowing for different expansion orders across space/energy. Our spatial discretisation is specifically designed to use less memory than competing DG schemes and also gives us direct access to the amount of stabilisation applied at each node. For filtered P〈sub〉〈em〉n〈/em〉〈/sub〉 expansions, we then use our adaptive process in combination with this net amount of stabilisation to compute a spatially dependent filter strength that does not depend on 〈em〉a priori〈/em〉 spatial information. This applies heavy filtering only where discontinuities are present, allowing the filtered P〈sub〉〈em〉n〈/em〉〈/sub〉 expansion to retain high-order convergence where possible. Regular and goal-based error metrics are shown and both the adapted P〈sub〉〈em〉n〈/em〉〈/sub〉 and adapted filtered P〈sub〉〈em〉n〈/em〉〈/sub〉 methods show significant reductions in DOFs and runtime. The adapted filtered P〈sub〉〈em〉n〈/em〉〈/sub〉 with our spatially dependent filter shows close to fixed iteration counts and up to high-order is even competitive with P〈sup〉0〈/sup〉 discretisations in problems with heavy advection.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Jingwei Hu, Shi Jin, Ruiwen Shu〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The Boltzmann equation may contain uncertainties in initial/boundary data or collision kernel. To study the impact of these uncertainties, a stochastic Galerkin (sG) method was proposed in [18] and studied in the kinetic regime. When the system is close to the fluid regime (the Knudsen number is small), the method would become prohibitively expensive due to the stiff collision term. In this work, we develop efficient sG methods for the Boltzmann equation that work for a wide range of Knudsen numbers, and investigate, in particular, their behavior in the fluid regime.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Qing Pan, Timon Rabczuk, Gang Xu, Chong Chen〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We investigate the isogeometric analysis for surface PDEs based on the extended Loop subdivision approach. The basis functions consisting of quartic box-splines corresponding to each subdivided control mesh are utilized to represent the geometry exactly, and construct the solution space for dependent variables as well, which is consistent with the concept of isogeometric analysis. The subdivision process is equivalent to the 〈em〉h〈/em〉-refinement of NURBS-based isogeometric analysis. The performance of the proposed method is evaluated by solving various surface PDEs, such as surface Laplace-Beltrami harmonic/biharmonic/triharmonic equations, which are defined on the limit surfaces of extended Loop subdivision for different initial control meshes. Numerical experiments show that the proposed method has desirable performance in terms of the accuracy, convergence and computational cost for solving the above surface PDEs defined on both open and closed surfaces. The proposed approach is proved to be second-order accuracy in the sense of 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈msup〉〈mrow〉〈mi〉L〈/mi〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈/mrow〉〈/msup〉〈/math〉-norm with theoretical and/or numerical results, which is also outperformed over the standard linear finite element by several numerical comparisons.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Arthur E.P. Veldman〈/p〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Massimiliano Ferronato, Andrea Franceschini, Carlo Janna, Nicola Castelletto, Hamdi A. Tchelepi〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This work discusses a general approach for preconditioning the block Jacobian matrix arising from the discretization and linearization of coupled multiphysics problem. The objective is to provide a fully algebraic framework that can be employed as a starting point for the development of specialized algorithms exploiting unique features of the specific problem at hand. The basic idea relies on approximately computing an operator able to decouple the different processes, which can then be solved independently one from the other. In this work, the decoupling operator is computed by extending the theory of block sparse approximate inverses. The proposed approach is implemented for two multiphysics applications, namely the simulation of a coupled poromechanical system and the mechanics of fractured media. The numerical results obtained in experiments taken from real-world examples are used to analyze and discuss the properties of the preconditioner.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 16 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Edoardo Zoni, Yaman Güçlü〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉A common strategy in the numerical solution of partial differential equations is to define a uniform discretization of a tensor-product multi-dimensional logical domain, which is mapped to a physical domain through a given coordinate transformation. By extending this concept to a multi-patch setting, simple and efficient numerical algorithms can be employed on relatively complex geometries. The main drawback of such an approach is the inherent difficulty in dealing with singularities of the coordinate transformation.〈/p〉 〈p〉This work suggests a comprehensive numerical strategy for the common situation of disk-like domains with a singularity at a unique pole, where one edge of the rectangular logical domain collapses to one point of the physical domain (for example, a circle). We present robust numerical methods for the solution of Vlasov-like hyperbolic equations coupled to Poisson-like elliptic equations in such geometries. We describe a semi-Lagrangian advection solver that employs a novel set of coordinates, named pseudo-Cartesian coordinates, to integrate the characteristic equations in the whole domain, including the pole, and a finite element elliptic solver based on globally 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈msup〉〈mrow〉〈mi mathvariant="script"〉C〈/mi〉〈/mrow〉〈mrow〉〈mn〉1〈/mn〉〈/mrow〉〈/msup〉〈/math〉 smooth splines (Toshniwal et al., 2017). The two solvers are tested both independently and on a coupled model, namely the 2D guiding-center model for magnetized plasmas, equivalent to a vorticity model for incompressible inviscid Euler fluids. The numerical methods presented show high-order convergence in the space discretization parameters, uniformly across the computational domain, without effects of order reduction due to the singularity. Dedicated tests show that the numerical techniques described can be applied straightforwardly also in the presence of point charges (equivalently, point-like vortices), within the context of particle-in-cell methods.〈/p〉 〈/div〉

Description: 〈p〉Publication date: Available online 16 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Maria Giuseppina Chiara Nestola, Barna Becsek, Hadi Zolfaghari, Patrick Zulian, Dario De Marinis, Rolf Krause, Dominik Obrist〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉We present a novel framework inspired by the Immersed Boundary Method for predicting the fluid-structure interaction of complex structures immersed in laminar, transitional and turbulent flows.〈/p〉 〈p〉The key elements of the proposed fluid-structure interaction framework are 1) the solution of elastodynamics equations for the structure, 2) the use of a high-order Navier–Stokes solver for the flow, and 3) the variational transfer (〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈msup〉〈mrow〉〈mi〉L〈/mi〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈/mrow〉〈/msup〉〈/math〉-projection) for coupling the solid and fluid subproblems.〈/p〉 〈p〉The dynamic behavior of a deformable structure is simulated in a finite element framework by adopting a fully implicit scheme for its temporal integration. It allows for mechanical constitutive laws including inhomogeneous and fiber-reinforced materials.〈/p〉 〈p〉The Navier–Stokes equations for the incompressible flow are discretized with high-order finite differences which allow for the direct numerical simulation of laminar, transitional and turbulent flows.〈/p〉 〈p〉The structure and the flow solvers are coupled by using an 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈msup〉〈mrow〉〈mi〉L〈/mi〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈/mrow〉〈/msup〉〈/math〉-projection method for the transfer of velocities and forces between the fluid grid and the solid mesh. This strategy allows for the numerical solution of coupled large scale problems based on nonconforming structured and unstructured grids. The transfer between fluid and solid limits the convergence order of the flow solver close to the fluid-solid interface.〈/p〉 〈p〉The framework is validated with the Turek–Hron benchmark and a newly proposed benchmark modelling the flow-induced oscillation of an inert plate. A three-dimensional simulation of an elastic beam in transitional flow is provided to show the solver's capability of coping with anisotropic elastic structures immersed in complex fluid flow.〈/p〉 〈/div〉

Description: 〈p〉Publication date: Available online 16 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Min Chai, Kun Luo, Changxiao Shao, Haiou Wang, Jianren Fan〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This paper proposes a finite difference discretization method for simulations of heat and mass transfer with Robin boundary conditions on irregular domains. The level set method is utilized to implicitly capture the irregular evolving interface, and the ghost fluid method to address variable discontinuities on the interface. Special care has been devoted to providing ghost values that are restricted by the Robin boundary conditions. Specifically, it is done in two steps: 1) calculate the normal derivative in cells adjacent to the interface by reconstructing a linear polynomial system; 2) successively extrapolate the normal derivative and the ghost value in the normal direction using a linear partial differential equation approach. This method produces second-order accurate solutions for both the Poisson and heat equations with Robin boundary conditions, and first-order accurate solutions for the Stefan problems. The solution gradients are of first-order accuracy, as expected. It is easy to implement in three-dimensional configurations, and can be straightforwardly generalized into higher-order variants. The method thus represents a promising tool for practical heat and mass transfer problems involving Robin boundary conditions.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Weidong Li, Wei Li, Pai Song, Hao Ji〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This paper presents an efficient, low memory cost, implicit finite volume lattice Boltzmann method (FVLBM) based on conservation moments acceleration for steady nearly incompressible flows. In the proposed scheme, not as the conventional implicit schemes, both the micro lattice Boltzmann equations (LBE) and the associated conservation moment equations are solved by the matrix-free, lower-upper symmetric Gauss-Seidel scheme (LUSGS) and the conservation moment equations are used to predict equilibrium distribution functions at the new time, which eliminates the storage of the Jacobian matrix of the collision term in the implicit LBE system and provides a driving force for the fast convergence of the LBE. Moreover, by utilizing the projection matrix and the collision invariant, we can construct the fluxes of the moment equations efficiently from the fluxes of the LBE and avoid the time-consuming reconstruction procedure for obtaining the fluxes of the moment equations. To demonstrate the accuracy and high efficiency of the proposed scheme, comparison studies of simulation results of several two-dimensional testing cases by the present scheme and an explicit FVLBM are conducted and numerical results show that the proposed implicit scheme can be as accurate as its explicit counterpart with 1 or 2 orders times speedup.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Laurent Muscat, Guillaume Puigt, Marc Montagnac, Pierre Brenner〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This paper addresses how two time integration schemes, the Heun's scheme for explicit time integration and the second-order Crank-Nicolson scheme for implicit time integration, can be coupled spatially. This coupling is the prerequisite to perform a coupled Large Eddy Simulation/Reynolds Averaged Navier-Stokes computation in an industrial context, using the implicit time procedure for the boundary layer (RANS) and the explicit time integration procedure in the LES region. The coupling procedure is designed in order to switch from explicit to implicit time integrations as fast as possible, while maintaining stability. After introducing the different schemes, the paper presents the initial coupling procedure adapted from a published reference and shows that it can amplify some numerical waves. An alternative procedure, studied in a coupled time/space framework, is shown to be stable and with spectral properties in agreement with the requirements of industrial applications. The coupling technique is validated with standard test cases, ranging from one-dimensional to three-dimensional flows.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Benjamin Stadlbauer, Andrea Cossettini, José A. Morales E., Daniel Pasterk, Paolo Scarbolo, Leila Taghizadeh, Clemens Heitzinger, Luca Selmi〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉Massively parallel nanosensor arrays fabricated with low-cost CMOS technology represent powerful platforms for biosensing in the Internet-of-Things (IoT) and Internet-of-Health (IoH) era. They can efficiently acquire “big data” sets of dependable calibrated measurements, representing a solid basis for statistical analysis and parameter estimation.〈/p〉 〈p〉In this paper we propose Bayesian estimation methods to extract physical parameters and interpret the statistical variability in the measured outputs of a dense nanocapacitor array biosensor. Firstly, the physical and mathematical models are presented. Then, a simple 1D-symmetry structure is used as a validation test case where the estimated parameters are also known 〈em〉a-priori〈/em〉. Finally, we apply the methodology to the simultaneous extraction of multiple physical and geometrical parameters from measurements on a CMOS pixelated nanocapacitor biosensor platform.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Meiling Zhao, Na Zhu〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We study the problem of electromagnetic scattering by multiple open cavities embedded in an infinite ground plane with high wave numbers. The problem can be described by a series of Helmholtz equations with coupled boundary conditions. We develop a sixth-order finite difference scheme to discretize the coupled Helmholtz equations. By Gaussian elimination in the vertical direction and Fourier transform in the horizontal direction, we can reduce the multiple cavity scattering problem to an aperture linear system. However, in the situation of high wave numbers, the condition number of the coefficient matrix of the reduced linear system is especially large and the system tends to be ill-conditioned. The convergence histories of most iterative methods become oscillating which consume considerable computations and memory spaces. In order to overcome the difficulty caused by high wave numbers, we develop an efficient preconditioned iterative method based on the Krylov subspace, which greatly improves the eigenvalue distributions and reduces the number of iterations. Numerical experiments show the validity and efficiency of the proposed sixth-order fast preconditioned algorithm for solving the scattering by multiple cavities with high wave numbers.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): YuFeng Shi, Yan Guo〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this paper, we propose an alternative formulation of conservative fifth order finite difference compact Weighted Essentially Non-Oscillatory (WENO) schemes to solve compressible Euler equations. Comparing with the classical conservative finite difference Compact-WENO scheme, its reconstruction procedure is applied to the point values rather than the traditional flux functions, then the HLLC (Harten, Lax and van Leer) and the local Lax-Friedrichs flux functions can be used to compute the interface fluxes in this framework. To maintain positivity of density and pressure, the parametrized positivity satisfying flux limiter is coupled with the proposed scheme for problems with extreme conditions. A number of testing cases including Titarev-Toro problem, the planar Sedov blast-wave problem, Riemann problems, double Mach reflection problem, shock diffraction problem and Kelvin-Helmholtz instability problem are presented to demonstrate the high resolution of the proposed compact scheme.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Alexandre Noll Marques, Jean-Christophe Nave, Rodolfo Ruben Rosales〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We introduce a technique that simplifies the problem of imposing jump conditions on interfaces that are not aligned with a computational grid in the context of the 〈em〉Correction Function Method〈/em〉 (CFM). The CFM offers a general framework to solve Poisson's equation in the presence of discontinuities to high order of accuracy, while using a compact discretization stencil. A key concept behind the CFM is enforcing the jump conditions in a least squares sense. This concept requires computing integrals over sections of the interface, which is a challenge in 3-D when only an implicit representation of the interface is available (e.g., the zero contour of a level set function). The technique introduced here is based on a new formulation of the least squares procedure that relies only on integrals over domains that are amenable to simple quadrature after local coordinate transformations. We incorporate this technique into a fourth order accurate implementation of the CFM, and show examples of solutions to Poisson's equation with imposed jump conditions computed in 2-D and 3-D.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Qian Kong, Yan-Fei Jing, Ting-Zhu Huang, Heng-Bin An〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The main aim of this paper is to develop two algorithms based on the Scheduled Relaxation Jacobi (SRJ) method (Yang and Mittal (2014) [7]) for solving problems arising from the finite-difference discretization of elliptic partial differential equations on large grids. These two algorithms are the Alternating Anderson-Scheduled Relaxation Jacobi (AASRJ) method by utilizing Anderson mixing after each SRJ iteration cycle and the Minimal Residual Scheduled Relaxation Jacobi (MRSRJ) method by minimizing residual after each SRJ iteration cycle, respectively. Through numerical experiments, we show that AASRJ is competitive with the optimal version of the SRJ method (Adsuara et al. (2017) [9]) in most problems we considered here, and MRSRJ outperforms SRJ in all cases. The properties of AASRJ and MRSRJ are demonstrated. Both of them are promising strategies for solving large, sparse linear systems while maintaining the simplicity of the Jacobi method.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Karim Shawki, George Papadakis〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We propose a preconditioner that can accelerate the rate of convergence of the Multiple Shooting Shadowing (MSS) method [1]. This recently proposed method can be used to compute derivatives of time-averaged objectives (also known as sensitivities) to system parameter(s) for chaotic systems. We propose a block diagonal preconditioner, which is based on a partial singular value decomposition of the MSS constraint matrix. The preconditioner can be computed using matrix-vector products only (i.e. it is matrix-free) and is fully parallelised in the time domain. Two chaotic systems are considered, the Lorenz system and the 1D Kuramoto Sivashinsky equation. Combination of the preconditioner with a regularisation method leads to tight bracketing of the eigenvalues to a narrow range. This combination results in a significant reduction in the number of iterations, and renders the convergence rate almost independent of the number of degrees of freedom of the system, and the length of the trajectory that is used to compute the time-averaged objective. This can potentially allow the method to be used for large chaotic systems (such as turbulent flows) and optimal control applications. The singular value decomposition of the constraint matrix can also be used to quantify the effect of the system condition on the accuracy of the sensitivities. In fact, neglecting the singular modes affected by noise, we recover the correct values of sensitivity that match closely with those obtained with finite differences for the Kuramoto Sivashinsky equation in the light turbulent regime. We notice a similar improvement for the Lorenz system as well.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Jian Yu, Chao Yan, Zhenhua Jiang, Wu Yuan, Shusheng Chen〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉An adaptive non-intrusive reduced basis (RB) method based on Gaussian process regression (GPR) is proposed for parametrized compressible flows. Adaptivity is pursued in the offline stage. The reduced basis by proper orthogonal decomposition (POD) is constructed iteratively to achieve a specified tolerance. For GPR, active data selection is used at each iteration, with standard deviation as the error indicator. To improve accuracy for shock-dominated flows, a properly designed simplified problem (SP) is considered as input of the regression models in addition to using parameters directly. Furthermore, a surrogate error model is constructed to serve as an efficient error estimator for the GPR models. Several two- and three-dimensional cases are conducted, including the inviscid nozzle flow, the inviscid NACA0012 airfoil flow and the inviscid M6 wing flow. For all the cases, the trained models are able to make efficient predictions with reasonable accuracy in the online stage. The SP-based approach is observed to result in biased sampling towards transonic regions. The regression models are further applied in sensitivity analysis, from which the solution of the two-dimensional cases are shown to be significantly more sensitive to input parameters than the wing flow. This is consistent to the comparison of convergence histories between the parameter-based and the SP-based models. For cases of high sensitivity, the SP-based approach is superior and can help to significantly reduce the number of required snapshots to achieve a prescribed tolerance.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): T. Allen, M. Zerroukat〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉With the increasing resolution of numerical weather prediction models the slopes of orographic features such as hills and mountains is becoming much better resolved with the consequence that the gradients encountered are becoming ever steeper. The commonly used terrain following coordinate, which actually becomes singular at 90°, is becoming more problematic due to the resulting stiffness making the elliptic boundary value problem difficult to solve. Even when a solution is obtainable the feedback of the pressure correction onto the momentum will generally cause instabilities and/or unphysical solutions. This paper uses an alternative representation of the orography, namely the immersed boundary method, in a semi-implicit-semi-Lagrangian model of the inviscid compressible Euler equations. The results show that the method not only performs well but, due to the chosen method of implementation, that the computational overhead of the geometry can be decoupled from implicit treatment of the vertical gravity wave component.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 31 July 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): François Forgues, Lucian Ivan, Alexandre Trottier, James G. M〈sup〉c〈/sup〉Donald〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The accurate prediction of multiphase flow when particles are differentiated by a set of “internal” variables, such as size or temperature, can pose modelling and numerical challenges. Although Lagrangian particle methods can provide predictive simulations of a wide spectrum of complex multiphase problems, they can become prohibitively expensive as the number of particles becomes large. Alternatively, Eulerian approaches have the potential to improve the computational efficiency of multiphase flows, but classical methods produce modelling artifacts or do not properly treat the local statistical dependence between the particle velocities and internal variables, or between the internal variables themselves. In this paper an extension is proposed of the classical Gaussian ten-moment model from gaskinetic theory to a model for the treatment of a dilute particle phase with an arbitrary number of internal variables based on an entropy-maximization argument. Unlike previous formulations, this new model provides a set of first-order robustly-hyperbolic balance laws that include a direct treatment for the local statistical variance of each variable, as well as the covariance between the internal variables or the internal variables and particle velocity. A study of the wave speeds of the general hyperbolic system is presented. To demonstrate an example application, the model is then specialized for polydisperse flows that are subject to viscous fluid drag as well as gravitational and buoyancy forces. The complete eigenstructure of this fifteen-moment polydisperse Gaussian model (PGM) is presented and the PGM is shown to maintain a physically realizable distribution function for all admissible initial conditions. Finally, several illustrative low-dimensional problems are studied to demonstrate the predictive capabilities of the new model.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 16 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): S. Danilov, A. Kutsenko〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉Computational wave branches are common to linearized shallow water equations discretized on triangular meshes. It is demonstrated that for standard finite-volume discretizations these branches can be traced back to the structure of the unit cell of triangular lattice, which includes two triangles with a common edge. Only subsets of similarly oriented triangles or edges possess the translational symmetry of unit cell. As a consequence, discrete degrees of freedom placed on triangles or edges are geometrically different, creating an internal structure inside unit cells. It implies a possibility of oscillations inside unit cells seen as computational branches in the framework of linearized shallow water equations, or as grid-scale noise generally.〈/p〉 〈p〉Adding dissipative operators based on smallest stencils to discretized equations is needed to control these oscillations in solutions. A review of several finite-volume discretization is presented with focus on computational branches and dissipative operators.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Claudio Canuto, Sandra Pieraccini, Dongbin Xiu〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A non-intrusive bifidelity strategy [1] is applied to the computation of statistics of a quantity of interest (QoI) which depends in a non-smooth way upon the stochastic parameters. The procedure leverages the accuracy of a high-fidelity model and the efficiency of a low-fidelity model, obtained through the use of different levels of numerical resolution, to pursue a high quality approximation of the statistics with a moderate number of high-fidelity simulations. The method is applied first to synthetic test cases with outputs exhibiting either a continuous or a discontinuous behaviour, then to the realistic simulation of a flow in an underground network of fractures [2], whose stochastic geometry outputs a non-smooth QoI. In both applications, the results highlight the efficacy of the approach in terms of error decay versus the number of computed high-fidelity solutions, even when the QoI lacks smoothness. For the underground simulation problem, the observed gain in computational cost is at least of one order of magnitude.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Samuel H. Rudy, Steven L. Brunton, J. Nathan Kutz〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The analysis of high-dimensional dynamical systems generally requires the integration of simulation data with experimental measurements. Experimental data often has substantial amounts of measurement noise that compromises the ability to produce accurate dimensionality reduction, parameter estimation, reduced order models, and/or balanced models for control. Data assimilation attempts to overcome the deleterious effects of noise by producing a set of algorithms for state estimation from noisy and possibly incomplete measurements. Indeed, methods such as Kalman filtering and smoothing are vital tools for scientists in fields ranging from electronics to weather forecasting. In this work we develop a novel framework for smoothing data based on known or partially known nonlinear governing equations. The method yields superior results to current techniques when applied to problems with known deterministic dynamics. By exploiting the numerical time-stepping constraints of the deterministic system, an optimization formulation can readily extract the noise from the nonlinear dynamics in a principled manner. The superior performance is due in part to the fact that it optimizes global state estimates. We demonstrate the efficiency and efficacy of the method on a number of canonical examples, thus demonstrating its viability for the wide range of potential applications stated above.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): M. Ganesh, S.C. Hawkins, D. Volkov〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We describe an efficient algorithm for reconstruction of the electromagnetic parameters of an unbounded dielectric medium from noisy cross section data induced by a point source in 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈msup〉〈mrow〉〈mi mathvariant="double-struck"〉R〈/mi〉〈/mrow〉〈mrow〉〈mn〉3〈/mn〉〈/mrow〉〈/msup〉〈/math〉. The efficiency of our Bayesian inverse algorithm for the parameters is based on developing an offline high order forward stochastic model and also an associated deterministic dielectric media Maxwell solver. Underlying the inverse/offline approach is our high order fully discrete Galerkin algorithm for solving an equivalent surface integral equation reformulation that is stable for all frequencies. The efficient algorithm includes approximating the likelihood distribution in the Bayesian model by a decomposed fast generalized polynomial chaos (gPC) model as a surrogate for the forward model. Offline construction of the gPC model facilitates fast online evaluation of the posterior distribution of the dielectric medium parameters. Parallel computational experiments demonstrate the efficiency of our deterministic, forward stochastic, and inverse dielectric computer models.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 8 August 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Jérôme Droniou, Matej Medla, Karol Mikula〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉We develop and analyse finite volume methods for the Poisson problem with boundary conditions involving oblique derivatives. We design a generic framework, for finite volume discretisations of such models, in which internal fluxes are not assumed to have a specific form, but only to satisfy some (usual) coercivity and consistency properties. The oblique boundary conditions are split into a normal component, which directly appears in the flux balance on control volumes touching the domain boundary, and a tangential component which is managed as an advection term on the boundary. This advection term is discretised using a finite volume method based on a centred discretisation (to ensure optimal rates of convergence) and stabilised using a vanishing boundary viscosity. A convergence analysis, based on the 3rd Strang Lemma [9], is conducted in this generic finite volume framework, and yields the expected 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈mi mathvariant="script"〉O〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈mi〉h〈/mi〉〈mo stretchy="false"〉)〈/mo〉〈/math〉 optimal convergence rate in discrete energy norm.〈/p〉 〈p〉We then describe a specific choice of numerical fluxes, based on a generalised hexahedral meshing of the computational domain. These fluxes are a corrected version of fluxes originally introduced in [29]. We identify mesh regularity parameters that ensure that these fluxes satisfy the required coercivity and consistency properties. The theoretical rates of convergence are illustrated by an extensive set of 3D numerical tests, including some conducted with two variants of the proposed scheme. A test involving real-world data measuring the disturbing potential in Earth gravity modelling over Slovakia is also presented.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Magnus Svärd, Jan Nordström〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉We consider constant-coefficient initial-boundary value problems, with a first or second derivative in time and a highest spatial derivative of order 〈em〉q〈/em〉, and their semi-discrete finite difference approximations. With an internal truncation error of order 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈mi〉p〈/mi〉〈mo〉≥〈/mo〉〈mn〉1〈/mn〉〈/math〉, and a boundary error of order 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.svg"〉〈mi〉r〈/mi〉〈mo〉≥〈/mo〉〈mn〉0〈/mn〉〈/math〉, we prove that the convergence rate is: 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si3.svg"〉〈mi mathvariant="normal"〉min〈/mi〉〈mo〉⁡〈/mo〉〈mo stretchy="false"〉(〈/mo〉〈mi〉p〈/mi〉〈mo〉,〈/mo〉〈mi〉r〈/mi〉〈mo linebreak="badbreak" linebreakstyle="after"〉+〈/mo〉〈mi〉q〈/mi〉〈mo stretchy="false"〉)〈/mo〉〈/math〉.〈/p〉 〈p〉The assumptions needed for these results to hold are: i) The continuous problem is linear and well-posed (with a smooth solution). ii) The numerical scheme is consistent, nullspace consistent, nullspace invariant, and energy stable. These assumptions are often satisfied for Summation-By-Parts schemes.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Vincent Mons, Qi Wang, Tamer A. Zaki〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Various ensemble-based variational (EnVar) data assimilation (DA) techniques are developed to reconstruct the spatial distribution of a scalar source in a turbulent channel flow resolved by direct numerical simulation (DNS). In order to decrease the computational cost of the DA procedure and improve its performance, Kriging-based interpolation is combined with EnVar DA, which enables the consideration of relatively large ensembles with moderate computational resources. The performance of the proposed Kriging-EnVar (KEnVar) DA scheme is assessed and favorably compared to that of standard EnVar and adjoint-based variational DA in various scenarios. Sparse regularization is implemented in the framework of EnVar DA in order to better tackle the case of concentrated scalar emissions. The problem of optimal sensor placement is also addressed, and it is shown that significant improvement in the quality of the reconstructed source can be obtained without supplementary computational cost once the ensemble required by the DA procedure is formed.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Jie Ding, Zhongming Wang, Shenggao Zhou〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉To study ion transport in electrolyte solutions, we propose numerical methods for a modified Poisson–Nernst–Planck model with ionic steric effects (SPNP). Positivity preserving schemes based on harmonic-mean approximations are proposed on a nonuniform mesh for the spatial discretization of the SPNP equations. Both explicit and semi-implicit discretization are considered in time. Numerical analysis shows that explicit forward Euler and semi-implicit trapezoidal discretization lead to schemes that maintain fully discrete solution positivity under a constraint on a mesh ratio, while the semi-implicit backward Euler discretization maintain fully discrete solution positivity 〈em〉unconditionally〈/em〉. We further study the condition numbers of the matrix associated with the semi-implicit backward Euler discretization, and establish an upper bound on condition numbers, indicating that the developed discretization based on harmonic-mean approximations effectively solves the issue of the presence of large condition numbers when using the Slotboom-type variables. Further numerical simulations confirm the analysis results on accuracy, positivity, and bounded condition numbers. The proposed schemes are also applied to study practical applications, such as the impact of self and cross steric interactions on ion distribution and rectifying behavior in a slit-shaped nanopore with surface charges. Possible extensions of the numerical methods to other modified PNP models with ion correlations and steric effects are also discussed.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Chinmay S. Kulkarni, Pierre F.J. Lermusiaux〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We propose a novel numerical methodology to compute the advective transport and diffusion-reaction of tracer quantities. The tracer advection occurs through flow map composition and is super-accurate, yielding numerical solutions almost devoid of compounding numerical errors, while allowing for direct parallelization in the temporal direction. It is computed by implicitly solving the characteristic evolution through a modified transport partial differential equation and domain decomposition in the temporal direction, followed by composition with the known initial condition. This advection scheme allows a rigorous computation of the spatial and temporal error bounds, yields an accuracy comparable to that of Lagrangian methods, and maintains the advantages of Eulerian schemes. We further show that there exists an optimal value of the composition timestep that yields the minimum total numerical error in the computations, and derive the expression for this value. We develop schemes for the addition of tracer diffusion, reaction, and source terms, and for the implementation of boundary conditions. Finally, the methodology is applied in three flow examples, namely an analytical reversible swirl flow, an idealized flow exiting a strait undergoing sudden expansion, and a realistic ocean flow in the Bismarck sea. New benchmark problems for advection-diffusion-reaction schemes are developed and used to compare and contrast results with those of classic schemes. The results highlight the theoretical properties of the methodology as well as its efficiency, super-accuracy with minimal numerical errors, and applicability in realistic simulations.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 December 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 398〈/p〉 〈p〉Author(s): Immo Huismann, Jörg Stiller, Jochen Fröhlich〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉High-order methods gain increased attention in computational fluid dynamics. However, due to the time step restrictions arising from the semi-implicit time stepping for the incompressible case, the potential advantage of these methods depends critically on efficient elliptic solvers. Due to the operation counts of operators scaling with the polynomial degree 〈em〉p〈/em〉 times the number of degrees of freedom 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈msub〉〈mrow〉〈mi〉n〈/mi〉〈/mrow〉〈mrow〉〈mi mathvariant="normal"〉DOF〈/mi〉〈/mrow〉〈/msub〉〈/math〉, the runtime of the best available multigrid solvers scales with 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.svg"〉〈mi mathvariant="script"〉O〈/mi〉〈mrow〉〈mo stretchy="true"〉(〈/mo〉〈mi〉p〈/mi〉〈mo〉⋅〈/mo〉〈msub〉〈mrow〉〈mi〉n〈/mi〉〈/mrow〉〈mrow〉〈mi mathvariant="normal"〉DOF〈/mi〉〈/mrow〉〈/msub〉〈mo stretchy="true"〉)〈/mo〉〈/mrow〉〈/math〉. This scaling with 〈em〉p〈/em〉 significantly lowers the applicability of high-order methods to high orders. While the operators for residual evaluation can be linearized when using static condensation, 〈span〉Schwarz〈/span〉-type smoothers require their inverses on fixed subdomains. No explicit inverse is known in the condensed case and matrix-matrix multiplications scale with 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si3.svg"〉〈mi〉p〈/mi〉〈mo〉⋅〈/mo〉〈msub〉〈mrow〉〈mi〉n〈/mi〉〈/mrow〉〈mrow〉〈mi mathvariant="normal"〉DOF〈/mi〉〈/mrow〉〈/msub〉〈/math〉. This paper derives a matrix-free explicit inverse for the static condensed operator in a cuboidal, Cartesian subdomain. It scales with 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si4.svg"〉〈msup〉〈mrow〉〈mi〉p〈/mi〉〈/mrow〉〈mrow〉〈mn〉3〈/mn〉〈/mrow〉〈/msup〉〈/math〉 per element, i.e. 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈msub〉〈mrow〉〈mi〉n〈/mi〉〈/mrow〉〈mrow〉〈mi mathvariant="normal"〉DOF〈/mi〉〈/mrow〉〈/msub〉〈/math〉 globally, and allows for a linearly scaling additive 〈span〉Schwarz〈/span〉 smoother, yielding a 〈em〉p〈/em〉-multigrid cycle with an operation count of 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si6.svg"〉〈mi mathvariant="script"〉O〈/mi〉〈mrow〉〈mo stretchy="true"〉(〈/mo〉〈msub〉〈mrow〉〈mi〉n〈/mi〉〈/mrow〉〈mrow〉〈mi mathvariant="normal"〉DOF〈/mi〉〈/mrow〉〈/msub〉〈mo stretchy="true"〉)〈/mo〉〈/mrow〉〈/math〉. The resulting solver uses fewer than four iterations for all polynomial degrees to reduce the residual by ten orders and has a runtime scaling linearly with 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈msub〉〈mrow〉〈mi〉n〈/mi〉〈/mrow〉〈mrow〉〈mi mathvariant="normal"〉DOF〈/mi〉〈/mrow〉〈/msub〉〈/math〉 for polynomial degrees at least up to 48. Furthermore the runtime is less than one microsecond per unknown over wide parameter ranges when using one core of a CPU, leading to time-stepping for the incompressible 〈span〉Navier-Stokes〈/span〉 equations using as much time for explicitly treated convection terms as for the elliptic solvers.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): J.E. Macías-Díaz, A.S. Hendy〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Departing from a two-dimensional hyperbolic system that describes the interaction between some activator and inhibitor substances in chemical reactions, we investigate a general form of that model using a finite-difference approach. The model under investigation is a nonlinear system consisting of two coupled partial differential equations with generalized reaction terms. The presence of two-dimensional diffusive terms consisting of fractional operators of the Riesz type is considered here, using spatial differentiation orders in the set 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.svg"〉〈mo stretchy="false"〉(〈/mo〉〈mn〉0〈/mn〉〈mo〉,〈/mo〉〈mn〉1〈/mn〉〈mo stretchy="false"〉)〈/mo〉〈mo〉∪〈/mo〉〈mo stretchy="false"〉(〈/mo〉〈mn〉1〈/mn〉〈mo〉,〈/mo〉〈mn〉2〈/mn〉〈mo stretchy="false"〉]〈/mo〉〈/math〉. We impose initial conditions on a closed and bounded rectangle, and a finite-difference methodology based on the use of fractional centered differences is proposed. Among the most important results of this work, we prove the existence and the uniqueness of the solutions of the numerical method, and establish analytically the second-order consistency of our scheme. Moreover, the discrete energy method is employed to prove the stability and the quadratic convergence of the technique. Some numerical simulations obtained through our method show the presence of Turing patterns and wave instabilities, in agreement with some reports found in the literature on superdiffusive hyperbolic activator-inhibitor systems. We show numerically that the presence of Turing patterns is independent of the size of the spatial domain. As a new application, we show that Turing patterns are also present in subdiffusive scenarios.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 November 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 397〈/p〉 〈p〉Author(s): Lisandro Dalcin, Diego Rojas, Stefano Zampini, David C. Del Rey Fernández, Mark H. Carpenter, Matteo Parsani〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present a novel technique for the imposition of non-linear entropy conservative and entropy stable solid wall boundary conditions for the compressible Navier–Stokes equations in the presence of an adiabatic wall, or a wall with a prescribed heat entropy flow. The procedure relies on the formalism and mimetic properties of diagonal-norm, summation-by-parts and simultaneous-approximation-term operators, and is a generalization of previous works on discontinuous interface coupling [1] and solid wall boundary conditions [2]. Using the method of lines, a semi-discrete entropy estimate for the entire domain is obtained when the proposed numerical imposition of boundary conditions are coupled with an entropy-conservative or entropy-stable discrete interior operator. The resulting estimate mimics the global entropy estimate obtained at the continuous level. The boundary data at the wall are weakly imposed using a penalty flux approach and a simultaneous-approximation-term technique for both the conservative variables and the gradient of the entropy variables. Discontinuous spectral collocation operators (mass lumped nodal discontinuous Galerkin operators), on high-order unstructured grids, are used for the purpose of demonstrating the robustness and efficacy of the new procedure for weakly enforcing boundary conditions. Numerical simulations confirm the non-linear stability of the proposed technique, with applications to three-dimensional subsonic and supersonic flows. The procedure described is compatible with any diagonal-norm summation-by-parts spatial operator, including finite element, finite difference, finite volume, discontinuous Galerkin, and flux reconstruction schemes.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): Takashi Shiroto, Naofumi Ohnishi, Yasuhiko Sentoku〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉For more than half a century, most of the plasma scientists have encountered a violation of the conservation laws of charge, momentum, and energy whenever they have numerically solved the first-principle equations of kinetic plasmas, such as the relativistic Vlasov–Maxwell system. This fatal problem is brought by the fact that both the Vlasov and Maxwell equations are indirectly associated with the conservation laws by means of some mathematical manipulations. Here we propose a quadratic conservative scheme, which can strictly maintain the conservation laws by discretizing the relativistic Vlasov–Maxwell system. A discrete product rule and summation-by-parts are the key players in the construction of the quadratic conservative scheme. Numerical experiments of the relativistic two-stream instability and relativistic Weibel instability prove the validity of our computational theory, and the proposed strategy will open the doors to the first-principle studies of mesoscopic and macroscopic plasma physics.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 14 August 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Mehdi Samiee, Mohsen Zayernouri, Mark M. Meerschaert〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present the stability and error analysis of the unified Petrov–Galerkin spectral method, developed in [1], for linear fractional partial differential equations with two-sided derivatives and constant coefficients in any (〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mn〉1〈/mn〉〈mo〉+〈/mo〉〈mi〉d〈/mi〉〈/math〉)-dimensional space-time hypercube, 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mi〉d〈/mi〉〈mo〉=〈/mo〉〈mn〉1〈/mn〉〈mo〉,〈/mo〉〈mn〉2〈/mn〉〈mo〉,〈/mo〉〈mn〉3〈/mn〉〈mo〉,〈/mo〉〈mo〉⋯〈/mo〉〈/math〉, subject to homogeneous Dirichlet initial/boundary conditions. Specifically, we prove the existence and uniqueness of the weak form and perform the corresponding stability and error analysis of the proposed method. Finally, we perform several numerical simulations to compare the theoretical and computational rates of convergence.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Giovanni Soligo, Alessio Roccon, Alfredo Soldati〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this work, we propose and test the validity of a modified Phase Field Method (PFM), which is specifically developed for large scale simulations of turbulent flows with large and deformable surfactant-laden droplets. The time evolution of the phase field, 〈em〉ϕ〈/em〉, and of the surfactant concentration field, 〈em〉ψ〈/em〉, are obtained from two Cahn–Hilliard-like equations together with a two-order-parameter Time-Dependent Ginzburg–Landau (TDGL) free energy functional. The modifications introduced circumvent existing limitations of current approaches based on PFM and improve the well-posedness of the model. The effect of surfactant on surface tension is modeled via an Equation Of State (EOS), further improving the flexibility of the approach. This method can efficiently handle topological changes, i.e. breakup and coalescence, and describe adsorption/desorption of surfactant. The capabilities of the proposed approach are tested in this paper against previous experimental results on the effects of surfactant on the deformation of a single droplet and on the interactions between two droplets. Finally, to appreciate the performances of the model on a large scale complex simulation, a qualitative analysis of the behavior of surfactant-laden droplets in a turbulent channel flow is presented and discussed.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Sergii V. Siryk〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We provide a careful Fourier analysis of the Guermond–Pasquetti mass lumping correction technique (Guermond and Pasquetti, 2013 [11]) applied to pure transport and convection–diffusion problems. In particular, it is found that increasing the number of corrections reduces the accuracy for problems with diffusion; however all the corrected schemes are more accurate than the consistent Galerkin formulation in this case. For the pure transport problems the situation is the opposite. We also investigate the differences between two numerical solutions – the consistent solution and the corrected ones, and show that increasing the number of corrections makes solutions of the corrected schemes closer to the consistent solution in all cases.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 377〈/p〉 〈p〉Author(s): Adam S. Jermyn〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Tensors are a natural way to express correlations among many physical variables, but storing tensors in a computer naively requires memory which scales exponentially in the rank of the tensor. This is not optimal, as the required memory is actually set not by the rank but by the mutual information amongst the variables in question. Representations such as the tensor tree perform near-optimally when the tree decomposition is chosen to reflect the correlation structure in question, but making such a choice is non-trivial and good heuristics remain highly context-specific. In this work I present two new algorithms for choosing efficient tree decompositions, independent of the physical context of the tensor. The first is a brute-force algorithm which produces optimal decompositions up to truncation error but is generally impractical for high-rank tensors, as the number of possible choices grows exponentially in rank. The second is a greedy algorithm, and while it is not optimal it performs extremely well in numerical experiments while having runtime which makes it practical even for tensors of very high rank.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 377〈/p〉 〈p〉Author(s): Lam H. Nguyen, Dominik Schillinger〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We describe a local iterative corrector scheme that significantly improves the accuracy of the multiscale finite element method (MsFEM). Our technique is based on the definition of a local corrector problem for each multiscale basis function that is driven by the residual of the previous multiscale solution. Each corrector problem results in a local corrector solution that improves the accuracy of the corresponding multiscale basis function at element interfaces. We cast the strategy of residual-driven correction in an iterative scheme that is straightforward to implement and, due to the locality of corrector problems, well-suited for parallel computing. We show that the iterative scheme converges to the best possible fine-mesh solution. Finally, we illustrate the effectiveness of our approach with multiscale benchmarks characterized by missing scale separation, including the microCT-based stress analysis of a vertebra with trabecular microstructure.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 377〈/p〉 〈p〉Author(s): Brody Bassett, Brian Kiedrowski〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The meshless local Petrov–Galerkin (MLPG) method is applied to the steady-state and 〈em〉k〈/em〉-eigenvalue neutron transport equations, which are discretized in energy using the multigroup approximation and in angle using the discrete ordinates approximation. To prevent oscillations in the neutron flux, the MLPG transport equation is stabilized by the streamline upwind Petrov–Galerkin (SUPG) method. Global neutron conservation is enforced by using moving least squares basis and weight functions and appropriate SUPG parameters. The cross sections in the transport equation are approximated in accordance with global particle balance and without constraint on their spatial dependence or the location of the basis and weight functions. The equations for the strong-form meshless collocation approach are derived for comparison to the MLPG equations. The method of manufactured solutions is used to verify the resulting MLPG method in one, two and three dimensions. Results for realistic problems, including two-dimensional pincells, a reflected ellipsoid and a three-dimensional problem with voids, are verified by comparison to Monte Carlo simulations.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 28 May 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Frederic Gibou, David Hyde, Ron Fedkiw〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present a review on numerical methods for simulating multiphase and free surface flows. We focus in particular on numerical methods that seek to preserve the discontinuous nature of the solutions across the interface between phases. We provide a discussion on the Ghost-Fluid and Voronoi Interface methods, on the treatment of surface tension forces that avoid stringent time step restrictions, on adaptive grid refinement techniques for improved efficiency and on parallel computing approaches. We present the results of some simulations obtained with these treatments in two and three spatial dimensions. We also provide a discussion of Machine Learning and Deep Learning techniques in the context of multiphase flows and propose several future potential research thrusts for using deep learning to enhance the study and simulation of multiphase flows.〈/p〉〈/div〉 〈h5〉Graphical abstract〈/h5〉 〈div〉〈p〉〈/p〉〈/div〉 〈figure〉〈img src="https://ars.els-cdn.com/content/image/1-s2.0-S0021999118303371-gr001.jpg" width="500" alt="Graphical abstract for this article" title=""〉〈/figure〉

Description: 〈p〉Publication date: Available online 26 October 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Xiaodong Liu, Jiguang Sun〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Inverse scattering has been an active research area for the past thirty years. While very successful in many cases, progress has lagged when only 〈em〉limited-aperture〈/em〉 measurement is available. In this paper, we perform some elementary study to recover data that can not be measured directly. In particular, we aim at recovering the 〈em〉full-aperture〈/em〉 far field data from 〈em〉limited-aperture〈/em〉 measurement. Due to the reciprocity relation, the multi-static response matrix (MSR) has a symmetric structure. Using the Green's formula and single layer potential, we propose two schemes to recover 〈em〉full-aperture〈/em〉 MSR. The recovered data is tested by a recently proposed direct sampling method and the factorization method. Numerical results show that it is possible to, at least, partially recover the missing data and consequently improve the reconstruction of the scatterer.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 377〈/p〉 〈p〉Author(s): D. Reiser, J. Romazanov, Ch. Linsmeier〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The problem of constructing Monte-Carlo solutions of drift-diffusion systems corresponding to Fokker–Planck equations with sources and sinks is revisited. Firstly, a compact formalism is introduced for the specific problem of stationary solutions. This leads to identification of the dwell time as the key quantity to characterize the system and to obtain a proper normalization for statistical analysis of numerical results. Secondly, the question of appropriate track length estimators for drift-diffusion systems is discussed for a 1D model system. It is found that a simple track length estimator can be given only for pure drift motion without diffusion. The stochastic nature of the diffusive part cannot be appropriately described by the path length of simulation particles. Further analysis of the usual situation with inhomogeneous drift and diffusion coefficients leads to an error estimate based on particle trajectories. The result for limits in grid cell size and time step used for the construction of Monte-Carlo trajectories resembles the Courant-Friedrichs-Lewy and von Neumann conditions for explicit methods.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Dinshaw S. Balsara, Roger Käppeli〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉The time-dependent equations of computational electrodynamics (CED) are evolved consistent with the divergence constraints on the electric displacement and magnetic induction vector fields. Respecting these constraints has proved to be very useful in the classic finite-difference time-domain (FDTD) schemes. As a result, there has been a recent effort to design finite volume time domain (FVTD) and discontinuous Galerkin time domain (DGTD) schemes that satisfy the same constraints and, nevertheless, draw on recent advances in higher order Godunov methods. This paper catalogues the first step in the design of globally constraint-preserving DGTD schemes. The algorithms presented here are based on a novel DG-like method that is applied to a Yee-type staggering of the electromagnetic field variables in the faces of the mesh. The other two novel building blocks of the method include constraint-preserving reconstruction of the electromagnetic fields and multidimensional Riemann solvers; both of which have been developed in recent years by the first author.〈/p〉 〈p〉The resulting DGTD scheme is linear, at least when limiters are not applied to the DG scheme. As a result, it is possible to carry out a von Neumann stability analysis of the entire suite of DGTD schemes for CED at orders of accuracy ranging from second to fourth. The analysis requires some simplifications in order to make it analytically tractable, however, it proves to be extremely instructive. A von Neumann stability analysis is a necessary precursor to the design of a full DGTD scheme for CED. It gives us the maximal CFL numbers that can be sustained by the DGTD schemes presented here at all orders. It also enables us to understand the wave propagation characteristics of the schemes in various directions on a Cartesian mesh. We find that constraint-preserving DGTD schemes permit CFL numbers that are competitive with conventional DG schemes. However, like conventional DG schemes, the CFL of DGTD schemes decreases with increasing order. To counteract that, we also present constraint-preserving PNPM schemes for CED. We find that the third and fourth order constraint-preserving DGTD and P1PM schemes have some extremely attractive properties when it comes to low-dispersion, low-dissipation propagation of electromagnetic waves in multidimensions. Numerical accuracy tests are also provided to support the von Neumann stability analysis. We expect these methods to play a role in those problems of engineering CED where exceptional precision must be achieved at any cost.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 15 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 377〈/p〉 〈p〉Author(s): Yinghe Qi, Jiacai Lu, Ruben Scardovelli, Stéphane Zaleski, Grétar Tryggvason〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In spite of considerable progress, computing curvature in Volume of Fluid (VOF) methods continues to be a challenge. The goal is to develop a function or a subroutine that returns the curvature in computational cells containing an interface separating two immiscible fluids, given the volume fraction in the cell and the adjacent cells. Currently, the most accurate approach is to fit a curve (2D), or a surface (3D), matching the volume fractions and finding the curvature by differentiation. Here, a different approach is examined. A synthetic data set, relating curvature to volume fractions, is generated using well-defined shapes where the curvature and volume fractions are easily found and then machine learning is used to fit the data (training). The resulting function is used to find the curvature for shapes not used for the training and implemented into a code to track moving interfaces. The results suggest that using machine learning to generate the relationship is a viable approach that results in reasonably accurate predictions.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Mani Razi, Robert M. Kirby, Akil Narayan〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this paper, we introduce a novel approach for the construction of multi-fidelity surrogate models with “discrete” fidelity levels. The notion of a discrete level of fidelity is in contrast to a mathematical model, for which the notion of refinement towards a high-fidelity model is relevant to sending a discretization parameter toward zero in a continuous way. Our notion of discrete fidelity levels encompasses cases for which there is no notion of convergence in terms of a fidelity parameter that can be sent to zero or infinity. The particular choice of how levels of fidelity are defined in this framework paves the way for using models that may have no apparent physical or mathematical relationship to the target high-fidelity model. However, our approach requires that models can produce results with a common set of parameters in the target model. Hence, fidelity level in this work is not directly representative of the degree of similarity of a low-fidelity model to a target high-fidelity model. In particular, we show that our approach is applicable to competitive ecological systems with different numbers of species, discrete-state Markov chains with a different number of states, polymer networks with a different number of connections, and nano-particle plasmonic arrays with a different number of scatterers. The results of this study demonstrate that our procedure boasts computational efficiency and accuracy for a wide variety of models and engineering systems.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Hasan Almanasreh〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this work we will treat the spurious eigenvalues obstacle that appears in the computation of the radial Dirac eigenvalue problem using numerical methods. The treatment of the spurious solution is based on applying Petrov–Galerkin finite element method. The significance of this work is the employment of just continuous basis functions, thus the need of a continuous function which has a continuous first derivative as a basis, as in [2], [3], is no longer required. The Petrov–Galerkin finite element method for the Dirac eigenvalue problem strongly depends on a stability parameter, 〈em〉τ〈/em〉, that controls the size of the diffusion terms added to the finite element formulation for the problem. The mesh-dependent parameter 〈em〉τ〈/em〉 is derived based on the given problem with the particular basis functions.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): S.B. Adrian, F.P. Andriulli, T.F. Eibert〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present a Calderón preconditioner for the electric field integral equation (EFIE), which does not require a barycentric refinement of the mesh and which yields a Hermitian, positive definite (HPD) system matrix allowing for the usage of the conjugate gradient (CG) solver. The resulting discrete equation system is immune to the low-frequency and the dense-discretization breakdown and, in contrast to existing Calderón preconditioners, no second discretization of the EFIE operator with Buffa–Christiansen (BC) functions is necessary. This preconditioner is obtained by leveraging on spectral equivalences between (scalar) integral operators, namely the single layer and the hypersingular operator known from electrostatics, on the one hand, and the Laplace–Beltrami operator on the other hand. Since our approach incorporates Helmholtz projectors, there is no search for global loops necessary and thus our method remains stable on multiply connected geometries. The numerical results demonstrate the effectiveness of this approach for both canonical and realistic (multi-scale) problems.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 21 February 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Mehdi Samiee, Mohsen Zayernouri, Mark M. Meerschaert〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We develop a unified Petrov–Galerkin spectral method for a class of fractional partial differential equations with two-sided derivatives and constant coefficients of the form 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mmultiscripts〉〈mrow〉〈mi mathvariant="script"〉D〈/mi〉〈/mrow〉〈mrow〉〈mi〉t〈/mi〉〈/mrow〉〈none〉〈/none〉〈none〉〈/none〉〈mrow〉〈mn〉2〈/mn〉〈mi〉τ〈/mi〉〈/mrow〉〈mprescripts〉〈/mprescripts〉〈mrow〉〈mn〉0〈/mn〉〈/mrow〉〈none〉〈/none〉〈/mmultiscripts〉〈mi〉u〈/mi〉〈mo〉+〈/mo〉〈msubsup〉〈mrow〉〈mo〉∑〈/mo〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈mo〉=〈/mo〉〈mn〉1〈/mn〉〈/mrow〉〈mrow〉〈mi〉d〈/mi〉〈/mrow〉〈/msubsup〉〈mo stretchy="false"〉[〈/mo〉〈msub〉〈mrow〉〈msub〉〈mrow〉〈mi〉c〈/mi〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉l〈/mi〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msub〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉a〈/mi〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msub〉〈msubsup〉〈mrow〉〈mi mathvariant="script"〉D〈/mi〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉x〈/mi〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈msub〉〈mrow〉〈mi〉μ〈/mi〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msubsup〉〈mi〉u〈/mi〉〈mo〉+〈/mo〉〈msub〉〈mrow〉〈msub〉〈mrow〉〈mi〉c〈/mi〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉r〈/mi〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msub〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉x〈/mi〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msub〉〈msubsup〉〈mrow〉〈mi mathvariant="script"〉D〈/mi〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉b〈/mi〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈msub〉〈mrow〉〈mi〉μ〈/mi〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msubsup〉〈mi〉u〈/mi〉〈mo stretchy="false"〉]〈/mo〉〈mo〉+〈/mo〉〈mi〉γ〈/mi〉〈mspace width="0.2em"〉〈/mspace〉〈mspace width="0.2em"〉〈/mspace〉〈mi〉u〈/mi〉〈mo〉=〈/mo〉〈msubsup〉〈mrow〉〈mo〉∑〈/mo〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈mo〉=〈/mo〉〈mn〉1〈/mn〉〈/mrow〉〈mrow〉〈mi〉d〈/mi〉〈/mrow〉〈/msubsup〉〈mo stretchy="false"〉[〈/mo〉〈msub〉〈mrow〉〈msub〉〈mrow〉〈mi〉κ〈/mi〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉l〈/mi〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msub〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉a〈/mi〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msub〉〈msubsup〉〈mrow〉〈mi mathvariant="script"〉D〈/mi〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉x〈/mi〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈msub〉〈mrow〉〈mi〉ν〈/mi〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msubsup〉〈mi〉u〈/mi〉〈mo〉+〈/mo〉〈msub〉〈mrow〉〈msub〉〈mrow〉〈mi〉κ〈/mi〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉r〈/mi〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msub〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉x〈/mi〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msub〉〈msubsup〉〈mrow〉〈mi mathvariant="script"〉D〈/mi〉〈/mrow〉〈mrow〉〈msub〉〈mrow〉〈mi〉b〈/mi〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈msub〉〈mrow〉〈mi〉ν〈/mi〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈/mrow〉〈/msub〉〈/mrow〉〈/msubsup〉〈mi〉u〈/mi〉〈mo stretchy="false"〉]〈/mo〉〈mo〉+〈/mo〉〈mi〉f〈/mi〉〈/math〉, where 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mn〉2〈/mn〉〈mi〉τ〈/mi〉〈mo〉∈〈/mo〉〈mo stretchy="false"〉(〈/mo〉〈mn〉0〈/mn〉〈mo〉,〈/mo〉〈mn〉2〈/mn〉〈mo stretchy="false"〉)〈/mo〉〈/math〉, 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si3.gif" overflow="scroll"〉〈mn〉2〈/mn〉〈mi〉τ〈/mi〉〈mo〉≠〈/mo〉〈mn〉1〈/mn〉〈/math〉, 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si4.gif" overflow="scroll"〉〈mn〉2〈/mn〉〈msub〉〈mrow〉〈mi〉μ〈/mi〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈/mrow〉〈/msub〉〈mo〉∈〈/mo〉〈mo stretchy="false"〉(〈/mo〉〈mn〉0〈/mn〉〈mo〉,〈/mo〉〈mn〉1〈/mn〉〈mo stretchy="false"〉)〈/mo〉〈/math〉 and 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si5.gif" overflow="scroll"〉〈mn〉2〈/mn〉〈msub〉〈mrow〉〈mi〉ν〈/mi〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈/mrow〉〈/msub〉〈mo〉∈〈/mo〉〈mo stretchy="false"〉(〈/mo〉〈mn〉1〈/mn〉〈mo〉,〈/mo〉〈mn〉2〈/mn〉〈mo stretchy="false"〉)〈/mo〉〈/math〉, in a (〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si6.gif" overflow="scroll"〉〈mn〉1〈/mn〉〈mo〉+〈/mo〉〈mi〉d〈/mi〉〈/math〉)-dimensional 〈em〉space–time〈/em〉 hypercube, 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si7.gif" overflow="scroll"〉〈mi〉d〈/mi〉〈mo〉=〈/mo〉〈mn〉1〈/mn〉〈mo〉,〈/mo〉〈mn〉2〈/mn〉〈mo〉,〈/mo〉〈mn〉3〈/mn〉〈mo〉,〈/mo〉〈mo〉⋯〈/mo〉〈/math〉, subject to homogeneous Dirichlet initial/boundary conditions. We employ the eigenfunctions of the fractional Sturm–Liouville eigen-problems of the first kind in [1], called 〈em〉Jacobi poly-fractonomial〈/em〉s, as temporal bases, and the eigen-functions of the boundary-value problem of the second kind as temporal test functions. Next, we construct our spatial basis/test functions using Legendre polynomials, yielding mass matrices being independent of the spatial fractional orders (〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si8.gif" overflow="scroll"〉〈msub〉〈mrow〉〈mi〉μ〈/mi〉〈/mrow〉〈mrow〉〈mi〉i〈/mi〉〈/mrow〉〈/msub〉〈mo〉,〈/mo〉〈mspace width="0.2em"〉〈/mspace〉〈msub〉〈mrow〉〈mi〉ν〈/mi〉〈/mrow〉〈mrow〉〈mi〉j〈/mi〉〈/mrow〉〈/msub〉〈mo〉,〈/mo〉〈mspace width="0.2em"〉〈/mspace〉〈mi〉i〈/mi〉〈mo〉,〈/mo〉〈mspace width="0.2em"〉〈/mspace〉〈mi〉j〈/mi〉〈mo〉=〈/mo〉〈mn〉1〈/mn〉〈mo〉,〈/mo〉〈mn〉2〈/mn〉〈mo〉,〈/mo〉〈mo〉⋯〈/mo〉〈mo〉,〈/mo〉〈mi〉d〈/mi〉〈/math〉). Furthermore, we formulate a novel unified fast linear solver for the resulting high-dimensional linear system based on the solution of generalized eigen-problem of spatial mass matrices with respect to the corresponding stiffness matrices, hence, making the complexity of the problem optimal, i.e., 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si9.gif" overflow="scroll"〉〈mi mathvariant="script"〉O〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈msup〉〈mrow〉〈mi〉N〈/mi〉〈/mrow〉〈mrow〉〈mi〉d〈/mi〉〈mo〉+〈/mo〉〈mn〉2〈/mn〉〈/mrow〉〈/msup〉〈mo stretchy="false"〉)〈/mo〉〈/math〉. We carry out several numerical test cases to examine the CPU time and convergence rate of the method. The corresponding stability and error analysis of the Petrov–Galerkin method are carried out in [2].〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Ruilian Du, Yubin Yan, Zongqi Liang〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A new high-order finite difference scheme to approximate the Caputo fractional derivative 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mfrac〉〈mrow〉〈mn〉1〈/mn〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈/mrow〉〈/mfrac〉〈mo stretchy="true" maxsize="2.4ex" minsize="2.4ex"〉(〈/mo〉〈mspace width="0.2em"〉〈/mspace〉〈mmultiscripts〉〈mrow〉〈mi〉D〈/mi〉〈/mrow〉〈mrow〉〈mi〉t〈/mi〉〈/mrow〉〈mrow〉〈mi〉α〈/mi〉〈/mrow〉〈mprescripts〉〈/mprescripts〉〈mrow〉〈mn〉0〈/mn〉〈/mrow〉〈mrow〉〈mi〉C〈/mi〉〈/mrow〉〈/mmultiscripts〉〈mi〉f〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈msub〉〈mrow〉〈mi〉t〈/mi〉〈/mrow〉〈mrow〉〈mi〉k〈/mi〉〈/mrow〉〈/msub〉〈mo stretchy="false"〉)〈/mo〉〈mo〉+〈/mo〉〈mspace width="0.2em"〉〈/mspace〉〈mmultiscripts〉〈mrow〉〈mi〉D〈/mi〉〈/mrow〉〈mrow〉〈mi〉t〈/mi〉〈/mrow〉〈mrow〉〈mi〉α〈/mi〉〈/mrow〉〈mprescripts〉〈/mprescripts〉〈mrow〉〈mn〉0〈/mn〉〈/mrow〉〈mrow〉〈mi〉C〈/mi〉〈/mrow〉〈/mmultiscripts〉〈mi〉f〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈msub〉〈mrow〉〈mi〉t〈/mi〉〈/mrow〉〈mrow〉〈mi〉k〈/mi〉〈mo〉−〈/mo〉〈mn〉1〈/mn〉〈/mrow〉〈/msub〉〈mo stretchy="false"〉)〈/mo〉〈mo stretchy="true" maxsize="2.4ex" minsize="2.4ex"〉)〈/mo〉〈/math〉, 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mi〉k〈/mi〉〈mo〉=〈/mo〉〈mn〉1〈/mn〉〈mo〉,〈/mo〉〈mn〉2〈/mn〉〈mo〉,〈/mo〉〈mo〉…〈/mo〉〈mo〉,〈/mo〉〈mi〉N〈/mi〉〈/math〉, with the convergence order 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si160.gif" overflow="scroll"〉〈mi〉O〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈mi mathvariant="normal"〉Δ〈/mi〉〈msup〉〈mrow〉〈mi〉t〈/mi〉〈/mrow〉〈mrow〉〈mn〉4〈/mn〉〈mo〉−〈/mo〉〈mi〉α〈/mi〉〈/mrow〉〈/msup〉〈mo stretchy="false"〉)〈/mo〉〈/math〉, 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si4.gif" overflow="scroll"〉〈mi〉α〈/mi〉〈mo〉∈〈/mo〉〈mo stretchy="false"〉(〈/mo〉〈mn〉1〈/mn〉〈mo〉,〈/mo〉〈mn〉2〈/mn〉〈mo stretchy="false"〉)〈/mo〉〈/math〉 is obtained when 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si123.gif" overflow="scroll"〉〈msup〉〈mrow〉〈mi〉f〈/mi〉〈/mrow〉〈mrow〉〈mo〉‴〈/mo〉〈/mrow〉〈/msup〉〈mo stretchy="false"〉(〈/mo〉〈msub〉〈mrow〉〈mi〉t〈/mi〉〈/mrow〉〈mrow〉〈mn〉0〈/mn〉〈/mrow〉〈/msub〉〈mo stretchy="false"〉)〈/mo〉〈mo〉=〈/mo〉〈mn〉0〈/mn〉〈/math〉, where Δ〈em〉t〈/em〉 denotes the time step size. Based on this scheme we introduce a finite difference method for solving fractional diffusion wave equation with the convergence order 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si6.gif" overflow="scroll"〉〈mi〉O〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈mi mathvariant="normal"〉Δ〈/mi〉〈msup〉〈mrow〉〈mi〉t〈/mi〉〈/mrow〉〈mrow〉〈mn〉4〈/mn〉〈mo〉−〈/mo〉〈mi〉α〈/mi〉〈/mrow〉〈/msup〉〈mo〉+〈/mo〉〈msup〉〈mrow〉〈mi〉h〈/mi〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈/mrow〉〈/msup〉〈mo stretchy="false"〉)〈/mo〉〈/math〉, where 〈em〉h〈/em〉 denotes the space step size. Numerical examples are given to show that the numerical results are consistent with the theoretical results.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 377〈/p〉 〈p〉Author(s): Jie Du, Yang Yang〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Local discontinuous Galerkin (LDG) methods are popular for convection-diffusion equations. In LDG methods, we introduce an auxiliary variable 〈em〉p〈/em〉 to represent the derivative of the primary variable 〈em〉u〈/em〉, and solve them on the same mesh. It is well known that the maximum-principle-preserving (MPP) LDG method is only available up to second-order accuracy. Recently, we introduced a new algorithm, and solve 〈em〉u〈/em〉 and 〈em〉p〈/em〉 on different meshes, and obtained stability and optimal error estimates. In this paper, we will continue this approach and construct MPP third-order LDG methods for convection-diffusion equations on overlapping meshes. The new algorithm is more flexible and does not increase any computational cost. Numerical evidence will be given to demonstrate the accuracy and good performance of the third-order MPP LDG method.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Yunchang Seol, Yu-Hau Tseng, Yongsam Kim, Ming-Chih Lai〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this paper, a two-dimensional immersed boundary method is developed to simulate the dynamics of Newtonian vesicle in viscoelastic Oldroyd-B fluid under shear flow. The viscoelasticity effect of extra stress is well incorporated into the immersed boundary formulation using the indicator function. Our numerical methodology is first validated in comparison with theoretical results in purely Newtonian fluid, and then a series of numerical experiments is conducted to study the effects of different dimensionless parameters on the vesicle motions. Although the tank-treading (TT) motion of Newtonian vesicle in Oldroyd-B fluid under shear flow can be observed just like in Newtonian fluid, it is surprising to find that the stationary inclination angle can be negative without the transition to tumbling (TB) motion. Moreover, the inertia effect plays a significant role that is able to turn the vesicle back to positive inclination angle through TT-TB-TT transition as the Reynolds number increases. To the best of our knowledge, this is the first numerical work for the detailed investigations of Newtonian vesicle dynamics suspended in viscoelastic Oldroyd-B fluid. We believe that our numerical results can be used to motivate further studies in theory and experiments for such coupling vesicle problems.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 378〈/p〉 〈p〉Author(s): Oscar P. Bruno, Martín Maas〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This paper introduces a fast algorithm, applicable throughout the electromagnetic spectrum, for the numerical solution of problems of scattering by periodic surfaces in two-dimensional space. The proposed algorithm remains highly accurate and efficient for challenging configurations including randomly rough surfaces, deep corrugations, large periods, near grazing incidences, and, importantly, Wood-anomaly resonant frequencies. The proposed approach is based on use of certain “shifted equivalent sources” which enable FFT acceleration of a Wood-anomaly-capable quasi-periodic Green function introduced recently (Bruno and Delourme (2014) [4]). The Green-function strategy additionally incorporates an exponentially convergent shifted version of the classical 〈em〉spectral〈/em〉 series for the Green function. While the computing-cost asymptotics depend on the asymptotic configuration assumed, the computing costs rise at most linearly with the size of the problem for a number of important rough-surface cases we consider. In practice, single-core runs in computing times ranging from a fraction of a second to a few seconds suffice for the proposed algorithm to produce highly-accurate solutions in some of the most challenging contexts arising in applications.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): Denis S. Grebenkov, Sergey D. Traytak〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The generalized method of separation of variables (GMSV) is applied to solve boundary value problems for the Laplace operator in three-dimensional domains with disconnected spherical boundaries (e.g., an arbitrary configuration of non-overlapping partially reactive spherical sinks or obstacles). We consider both exterior and interior problems and all most common boundary conditions: Dirichlet, Neumann, Robin, and conjugate one. Using the translational addition theorems for solid harmonics to switch between the local spherical coordinates, we obtain a semi-analytical expression of the Green function as a linear combination of partial solutions whose coefficients are fixed by boundary conditions. Although the numerical computation of the coefficients involves series truncation and solution of a system of linear algebraic equations, the use of the solid harmonics as basis functions naturally adapted to the intrinsic symmetries of the problem makes the GMSV particularly efficient, especially for exterior problems. The obtained Green function is the key ingredient to solve boundary value problems and to determine various characteristics of stationary diffusion such as reaction rate, escape probability, harmonic measure, residence time, and mean first passage time, to name but a few. The relevant aspects of the numerical implementation and potential applications in chemical physics, heat transfer, electrostatics, and hydrodynamics are discussed.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 377〈/p〉 〈p〉Author(s): Samar Chehade, Audrey Kamta Djakou, Michel Darmon, Gilles Lebeau〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Non Destructive Examination (NDE) of industrial structures requires the modeling of specimen geometry echoes generated by the surfaces (entry, backwall …) of inspected blocks. For that purpose, the study of plane wave diffraction by a wedge is of great interest. The work presented here is preliminary research to model the case of an elastic wave diffracted by a wedge in the future, for which there exist various modeling approaches but the numerical aspects have only been developed for wedge angles lower than 〈em〉π〈/em〉. The spectral functions method has previously been introduced to solve the 2D diffraction problem of an immersed elastic wedge for angles lower than 〈em〉π〈/em〉. As a first step, the spectral functions method has been developed here for the diffraction on an acoustic wave by a stress-free wedge, in 2D and for any wedge angle, before studying the elastic wave diffraction from a wedge. In this method, the solution to the diffraction problem is expressed in terms of two unknown functions called the spectral functions. These functions are computed semi-analytically, meaning that they are the sum of two terms. One of them is determined exactly and the other is approached numerically, using a collocation method. A successful numerical validation of the method for all wedge angles is proposed, by comparison with the GTD (Geometrical Theory of Diffraction) solution derived from the exact Sommerfeld integral.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 377〈/p〉 〈p〉Author(s): Hong Fang, Yikun Hu, Caihui Yu, Ming Tie, Jie Liu, Chunye Gong〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The mesh deformation method based on radial basis functions (RBF) has many advantages and is widely used. RBF based mesh deformation method mainly has two steps: data reduction and displacement interpolation. The data reduction step includes solving interpolation weight coefficients and searching for the node with the maximum interpolation error. The data reduction schemes based on greedy algorithm is used to select an optimum reduced set of surface mesh nodes. In this paper, a parallel mesh deformation method based on parallel data reduction and displacement interpolation is proposed. The proposed recurrence Choleskey decomposition method (RCDM) can decrease the computational cost of solving interpolation weight coefficients from 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mi〉O〈/mi〉〈mrow〉〈mo stretchy="true"〉(〈/mo〉〈msubsup〉〈mrow〉〈mi〉N〈/mi〉〈/mrow〉〈mrow〉〈mi〉c〈/mi〉〈/mrow〉〈mrow〉〈mn〉4〈/mn〉〈/mrow〉〈/msubsup〉〈mo stretchy="true"〉)〈/mo〉〈/mrow〉〈/math〉 to 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mi〉O〈/mi〉〈mrow〉〈mo stretchy="true"〉(〈/mo〉〈msubsup〉〈mrow〉〈mi〉N〈/mi〉〈/mrow〉〈mrow〉〈mi〉c〈/mi〉〈/mrow〉〈mrow〉〈mn〉3〈/mn〉〈/mrow〉〈/msubsup〉〈mo stretchy="true"〉)〈/mo〉〈/mrow〉〈/math〉, where 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si3.gif" overflow="scroll"〉〈msub〉〈mrow〉〈mi〉N〈/mi〉〈/mrow〉〈mrow〉〈mi〉c〈/mi〉〈/mrow〉〈/msub〉〈/math〉 denotes the number of support nodes. The technology of parallel computing is used to accelerate the searching for the node with the maximum interpolation error and displacement interpolation. The combination of parallel data reduction and parallel interpolation can greatly improve the efficiency of mesh deformation. Two typical deformation problems of the ONERA M6 and DLR-F6 wing-body-Nacelle-Pylon configuration are taken as the test cases to validate the proposed approach and can get up to 19.57 times performance improvement with the proposed approach. Finally, the aeroelastic response of HIRENASD wing-body configuration is used to verify the efficiency and robustness of the proposed method.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Daniil Bochkov, Frederic Gibou〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present two finite volume schemes to solve a class of Poisson-type equations subject to Robin boundary conditions in irregular domains with 〈em〉piecewise smooth〈/em〉 boundaries. The first scheme results in a symmetric linear system and produces second-order accurate numerical solutions with first-order accurate gradients in the 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈msup〉〈mrow〉〈mi〉L〈/mi〉〈/mrow〉〈mrow〉〈mo〉∞〈/mo〉〈/mrow〉〈/msup〉〈/math〉-norm (for solutions with two bounded derivatives). The second scheme is nonsymmetric but produces second-order accurate numerical solutions as well as second-order accurate gradients in the 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈msup〉〈mrow〉〈mi〉L〈/mi〉〈/mrow〉〈mrow〉〈mo〉∞〈/mo〉〈/mrow〉〈/msup〉〈/math〉-norm (for solutions with three bounded derivatives). Numerical examples are given in two and three spatial dimensions.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 19 October 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Mandeep Deka, Shuvayan Brahmachary, Ramakrishnan Thirumalaisamy, Amaresh Dalal, Ganesh Natarajan〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We describe a new and simple strategy based on the Gauss divergence theorem for obtaining centroidal gradients on unstructured meshes. Unlike the standard Green–Gauss (SGG) reconstruction which requires face values of quantities whose gradients are sought, the proposed approach reconstructs the gradients using the normal derivative(s) at the faces. The new strategy, referred to as the Modified Green–Gauss (MGG) reconstruction results in consistent gradients which are at least first-order accurate on arbitrary polygonal meshes. We show that the MGG reconstruction is linearity preserving independent of the mesh topology and retains the consistent behaviour of gradients even on meshes with large curvature and high aspect ratios. The gradient accuracy in MGG reconstruction depends on the accuracy of discretisation of the normal derivatives at faces and this necessitates an iterative approach for gradient computation on non-orthogonal meshes. Numerical studies on different mesh topologies demonstrate that MGG reconstruction gives accurate and consistent gradients on non-orthogonal meshes, with the number of iterations proportional to the extent of non-orthogonality. The MGG reconstruction is found to be consistent even on meshes with large aspect ratio and curvature with the errors being lesser than those from linear least-squares reconstruction. A non-iterative strategy in conjunction with MGG reconstruction is proposed for gradient computations in finite volume simulations that achieves the accuracy and robustness of MGG reconstruction at a cost equivalent to that of SGG reconstruction. The efficacy of this strategy for fluid flow problems is demonstrated through numerical investigations in both incompressible and compressible regimes. The MGG reconstruction may, therefore, be viewed as a novel and promising blend of least-squares and Green–Gauss based approaches which can be implemented with little effort in open-source finite-volume solvers and legacy codes.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 18 July 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Marvin Bohm, Andrew R. Winters, Gregor J. Gassner, Dominik Derigs, Florian Hindenlang, Joachim Saur〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉The first paper of this series presents a discretely entropy stable discontinuous Galerkin (DG) method for the resistive magnetohydrodynamics (MHD) equations on three-dimensional curvilinear unstructured hexahedral meshes. Compared to other fluid dynamics systems such as the shallow water equations or the compressible Navier–Stokes equations, the resistive MHD equations need special considerations because of the divergence-free constraint on the magnetic field. For instance, it is well known that for the symmetrization of the ideal MHD system as well as the continuous entropy analysis a non-conservative term proportional to the divergence of the magnetic field, typically referred to as the Powell term, must be included. As a consequence, the mimicry of the continuous entropy analysis in the discrete sense demands a suitable DG approximation of the non-conservative terms in addition to the ideal MHD terms.〈/p〉 〈p〉This paper focuses on the 〈em〉resistive〈/em〉 MHD equations: Our first contribution is a proof that the resistive terms are symmetric and positive-definite when formulated in entropy space as gradients of the entropy variables, which enables us to show that the entropy inequality holds for the resistive MHD equations. This continuous analysis is the key for our DG discretization and guides the path for the construction of an approximation that discretely mimics the entropy inequality, typically termed 〈em〉entropy stability〈/em〉. Our second contribution is a detailed derivation and analysis of the discretization on three-dimensional curvilinear meshes. The discrete analysis relies on the summation-by-parts property, which is satisfied by the DG spectral element method (DGSEM) with Legendre–Gauss–Lobatto (LGL) nodes. Although the divergence-free constraint is included in the non-conservative terms, the resulting method has no particular treatment of the magnetic field divergence errors, which might pollute the solution quality. Our final contribution is the extension of the standard resistive MHD equations and our DG approximation with a divergence cleaning mechanism that is based on a generalized Lagrange multiplier (GLM).〈/p〉 〈p〉As a conclusion to the first part of this series, we provide detailed numerical validations of our DGSEM method that underline our theoretical derivations. In addition, we show a numerical example where the entropy stable DGSEM demonstrates increased robustness compared to the standard DGSEM.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): James F. Kelly, Harish Sankaranarayanan, Mark M. Meerschaert〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This paper develops appropriate boundary conditions for the two-sided fractional diffusion equation, where the usual second derivative in space is replaced by a weighted average of positive (left) and negative (right) fractional derivatives. Mass preserving, reflecting boundary conditions for two-sided fractional diffusion involve a balance of left and right fractional derivatives at the boundary. Stable, consistent explicit and implicit Euler methods are detailed, and steady state solutions are derived. Steady state solutions for two-sided fractional diffusion equations using both Riemann–Liouville and Caputo flux are computed. For Riemann–Liouville flux and reflecting boundary conditions, the steady-state solution is singular at one or both of the end-points. For Caputo flux and reflecting boundary conditions, the steady-state solution is a constant function. Numerical experiments illustrate the convergence of these numerical methods. Finally, the influence of the reflecting boundary on the steady-state behavior subject to both the Riemann–Liouville and Caputo fluxes is discussed.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): N. Saini, C. Kleinstreuer〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉All natural and a growing number of manufactured solid particles are non-spherical. Interesting fluid–particle dynamics applications include the transport of granular material, piling of seeds or grains, inhalation of toxic aerosols, use of nanofluids for enhanced cooling or improved lubrication, and optimal drug-targeting of tumors. A popular approach for computer simulations of such scenarios is the multi-sphere (MS) method, where any non-spherical particle is represented by an assemblage of spheres. However, the MS approach may lead to multiple sphere-to-sphere contact points during collision, and subsequently to erroneous particle transport and deposition. In cases where non-spherical particles can be approximated as ellipsoids with arbitrary aspect ratios, a new theory for particle transport, collision and wall interaction is presented which is more accurate computationally and more efficient than the MS method. In general, with the new ellipsoidal particle interaction (EPI) model, contact points and planes of ellipsoids, rather than spheres, are obtained based on a geometric potential algorithm. Then, interaction forces and torques of the colliding particles are determined via inscribed ‘pseudo-spheres’, employing the soft-particle approach. The off-center forces and moments are then transferred to the mass center of the ellipsoids to solve the appropriate translatory and angular equations of motion. Considering ellipses to illustrate the workings and predictive power of the new collision model, turbulent fluid–particle flow with the EPI model in a 2-D channel is simulated and compared with 3-D numerical benchmark results which relied on the MS method. The 2-D concentrations of micron particles with different aspect ratios matched closely with the 3-D cases. However, interesting differences occurred when comparing the particle-velocity profiles for which the 2-D EPI model generated somewhat larger particle velocities due to out-of-plane collisions, slightly higher particle–wall interactions, and two-way coupling effects.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Kyongmin Yeo, Igor Melnyk〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present a deep learning model, DE-LSTM, for the simulation of a stochastic process with an underlying nonlinear dynamics. The deep learning model aims to approximate the probability density function of a stochastic process via numerical discretization and the underlying nonlinear dynamics is modeled by the Long Short-Term Memory (LSTM) network. It is shown that, when the numerical discretization is used, the function estimation problem can be solved by a multi-label classification problem. A penalized maximum log likelihood method is proposed to impose a smoothness condition in the prediction of the probability distribution. We show that the time evolution of the probability distribution can be computed by a high-dimensional integration of the transition probability of the LSTM internal states. A Monte Carlo algorithm to approximate the high-dimensional integration is outlined. The behavior of DE-LSTM is thoroughly investigated by using the Ornstein–Uhlenbeck process and noisy observations of nonlinear dynamical systems; Mackey–Glass time series and forced Van der Pol oscillator. It is shown that DE-LSTM makes a good prediction of the probability distribution without assuming any distributional properties of the stochastic process. For a multiple-step forecast of the Mackey–Glass time series, the prediction uncertainty, denoted by the 95% confidence interval, first grows, then dynamically adjusts following the evolution of the system, while in the simulation of the forced Van der Pol oscillator, the prediction uncertainty does not grow in time even for a 3,000-step forecast.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): Christopher Eldred, Thomas Dubos, Evaggelos Kritsikis〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The rotating shallow water (RSW) equations are the usual testbed for the development of numerical methods for three-dimensional atmospheric and oceanic models. However, an arguably more useful set of equations are the thermal shallow water equations (TSW), which introduce an additional thermodynamic scalar but retain the single layer, two-dimensional structure of the RSW. As a stepping stone towards a three-dimensional atmospheric dynamical core, this work presents a quasi-Hamiltonian discretization of the thermal shallow water equations using compatible Galerkin methods, building on previous work done for the shallow water equations. Structure-preserving or quasi-Hamiltonian discretizations methods, that discretize the Hamiltonian structure of the equations of motion rather than the equations of motion themselves, have proven to be a powerful tool for the development of models with discrete conservation properties. By combining these ideas with an energy-conserving Poisson time integrator and a careful choice of Galerkin spaces, a large set of desirable properties can be achieved. In particular, for the first time total mass, buoyancy and energy are conserved to machine precision in the fully discrete model.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 377〈/p〉 〈p〉Author(s): Fang Qing, Xijun Yu, Zupeng Jia〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The MoF (Moment of Fluid) method is an accurate approach for interface reconstruction in numerical simulation of multi-material fluid flow. So far, most works focus on improving its accuracy and efficiency, such as developing analytic reconstruction method and deducing the iteration schemes based on high order derivatives of the objective function. In this paper, we mainly concern on improving its robustness, especially for severely deformed polygonal meshes, in which case the objective function has multiple minimum value points. By using an efficient method for solving multiple roots of the nonlinear equation in large scope, a new algorithm is developed to enhance robustness of the MoF method. The main idea of this algorithm is as follows. The first derivative of the objective function is continuous, so the minimum value points of the objective function must be the zero points of the first derivative. Instead of finding the zero points of the first derivative directly, we turn to calculating the minimum value points (also zero points) of the square of the first derivative, which is a convex function on a neighborhood of each zero point. Applying the properties of convex function, the neighbor of each extreme minimum point of it can be obtained efficiently. Then each zero point of the square of the first derivative can be obtained using the iterative formula in its neighbor. Finally, by comparing the values of the objective function at these zero points of the first derivative, the global minimum value point of the objective function can be found and is the desired solution. The new algorithm only uses the first derivative of the objective function. It doesn't need an initial guess for the solution, which has to be carefully chosen in previous works. Numerical results are presented to demonstrate the accuracy and robustness of this new algorithm. The results show that it is applicable to severely deformed polygonal mesh, even with concave cells.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 22 March 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Savio Poovathingal, Eric C. Stern, Ioannis Nompelis, Thomas E. Schwartzentruber, Graham V. Candler〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Micro scale simulations are performed of flow through porous (pyrolyzing) thermal protection system (TPS) materials using the direct simulation Monte Carlo (DSMC) method. DSMC results for permeability are validated with computational fluid dynamics (CFD) calculations and theory, for simple porous geometries under continuum flow conditions. An artificial fiber-microstructure generation code FiberGen is used to create triangulated surface geometry representative of FiberForm® (FiberForm) material. DSMC results for permeability of FiberForm are validated for a range of pressures (transitional flow conditions) and agree with experimental measurements. Numerical uncertainty is determined to be within 2% if sufficiently large portions of the microstructure are included in the computation. However, small variations in fiber size and angle bias can combine to give 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mo〉+〈/mo〉〈mn〉30〈/mn〉〈mtext〉%〈/mtext〉〈/math〉 uncertainty when comparing with experimental permeability data. X-ray microtomography scans of FiberForm are used to create microstructure geometry for incorporation within DSMC simulations of coupled oxygen diffusion and gas-surface chemistry in the presence of a blowing pyrolysis gas. In-depth penetration of atomic oxygen is limited to 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mn〉0.2〈/mn〉〈mtext〉–〈/mtext〉〈mn〉0.4〈/mn〉〈/math〉 mm for the range of Knudsen number and pyrolysis gas conditions studied.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): A. Belme, F. Alauzet, A. Dervieux〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present a goal-oriented error analysis for the calculation of low Reynolds steady compressible flows with anisotropic mesh adaptation. The error analysis is of 〈em〉a priori〈/em〉 type. Its central principle is to express the right-hand side of the error equation, often referred as the local error, as a function of the interpolation error of a collection of fields present in the nonlinear Partial Differential Equations. This goal-oriented error analysis is the extension of [39] done for inviscid flows to laminar viscous flows by adding viscous terms. The main benefits of this approach, in comparison to other error estimates in the literature, is that the optimal anisotropy of the mesh directly appears in the error analysis and is not obtained from an ad hoc variable nor a local analysis. As a consequence, an optimum is obtained and the convergence of the mesh adaptive process is very fast, 〈em〉i.e.,〈/em〉 generally the convergence is obtained after 5 to 10 mesh adaptation cycle. Then, using the continuous mesh framework, an optimal metric is analytically obtained from the error estimation. Applications to mesh adaptive calculations of flows past airfoils are presented.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Pablo Miguel Ramos, Nikos Ch. Karayiannis, Manuel Laso〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A new algorithmic approach is presented for the generation and successive equilibration of polymer configurations under conditions of extreme confinement where the inter-wall distance, in at least one dimension, approaches the diameter of the spherical monomers. It significantly improves on the Monte Carlo (MC) protocol described in Karayiannis and Laso (2008) [126]. The algorithm is designed to generate highly confined packings of freely-jointed chains of hard spheres of uniform size. Spatial confinement is achieved by including flat, parallel impenetrable walls in one or more dimensions of the simulation box. The present MC scheme allows the systematic study of the effect of chain length, polydispersity, volume fraction, bond tolerance (gap), cell aspect ratio and level of confinement on the short- and long-range structure of polymer chains near and far from the confining planes. In the present study we focus on the efficiency of the MC protocol in generating, equilibrating, and configurationally decorrelating chain assemblies with average lengths ranging from 〈em〉N〈/em〉 = 12 to 1000 monomers and at volume fractions from dilute up to the maximally random jammed (MRJ) state. Starting from cubic amorphous cells filled with polymer chains, the MC algorithm is able to reach quasi 2-d (plate-like) and 1-d (tube-like) states under conditions of extreme confinement and/or cell aspect ratio where the inter-wall distance approaches the diameter of beads forming the chains. A comparison with corresponding bulk packings shows the similarities and differences produced by extreme spatial confinement.〈/p〉〈/div〉 〈h5〉Graphical abstract〈/h5〉 〈div〉〈p〉Snapshots of computer-generated athermal polymer configurations under full confinement with increasing cell aspect ratio starting from 3-d cubic cells and leading to 2-d templates.〈/p〉〈/div〉 〈figure〉〈img src="https://ars.els-cdn.com/content/image/1-s2.0-S0021999118305850-gr001.jpg" width="452" alt="Graphical abstract for this article" title=""〉〈/figure〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Gaddiel Y. Ouaknin, Nabil Laachi, Kris Delaney, Glenn H. Fredrickson, Frederic Gibou〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We introduce a level-set strategy to find the geometry of confinement that will guide the self-assembly of block copolymers to a given target design in the context of lithography. The methodology is based on a shape optimization algorithm, where the level-set normal velocity is defined as the pressure field computed through a self-consistent field theory simulation. We present numerical simulations that demonstrate that this methodology is capable of finding guiding templates for a variety of target arrangements of cylinders and thus is an effective approach to the inverse directed self-assembly problem.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Gautier Dakin, Bruno Després, Stéphane Jaouen〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Staggered grids schemes, formulated in internal energy, are commonly used for CFD applications in industrial context. Here, we prove the consistency of a class of high-order Lagrange-Remap staggered schemes for solving the Euler equations in 1D and 2D on Cartesian grids. The main result of the paper is that using an 〈em〉a posteriori〈/em〉 internal energy corrector, the Lagrangian schemes are proved to be conservative in mass, momentum and total energy and to be weakly consistent with the 1D Lagrangian formulation of the Euler equations. Extension in 2D is done using directional splitting methods and face-staggering. Numerical examples in both 1D and 2D illustrate the accuracy, the convergence and the robustness of the schemes.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Florian Monteghetti, Denis Matignon, Estelle Piot〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Time-domain impedance boundary conditions (TDIBCs) can be enforced using the impedance, the admittance, or the scattering operator. This article demonstrates the computational advantage of the last, even for nonlinear TDIBCs, with the linearized Euler equations. This is achieved by a systematic semi-discrete energy analysis of the weak enforcement of a generic nonlinear TDIBC in a discontinuous Galerkin finite element method. In particular, the analysis highlights that the sole definition of a discrete model is not enough to fully define a TDIBC. To support the analysis, an elementary physical nonlinear scattering operator is derived and its computational properties are investigated in an impedance tube. Then, the derivation of time-delayed broadband TDIBCs from physical reflection coefficient models is carried out for single degree of freedom acoustical liners. A high-order discretization of the derived time-local formulation, which consists in composing a set of ordinary differential equations with a transport equation, is applied to two flow ducts.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): C. Cheng, A.P. Bunger〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A new reduced order model (ROM) provides rapid and reasonably accurate prediction of the complex behavior of multiple, simultaneously growing radial hydraulic fractures. The method entails vastly reducing the degrees of freedom typically associated with fully-coupled simulations of this multiple moving boundary problem by coupling together an approximation of the influence of the stress interaction among the fractures (“stress shadow”) with an approximation of fluid flow and elasticity, ensuring preservation of global volume balance, global energy balance, elasticity, and compatibility of the crack opening with the inlet fluid flux. Validating with large scale (“high-fidelity”) simulations shows the ROM solution captures not only the basic suppression of interior hydraulic fractures in a uniformly-spaced array due to the well-known stress shadowing phenomenon, but also complex behaviors arising when the spacing among the hydraulic fractures is non-uniform. The simulator's usefulness is demonstrated through a proof-of-concept optimization whereby non-uniform spacing and stage length are chosen to maximize the fracture surface area and/or the uniformity of growth associated with each stimulation treatment.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Rohit K. Tripathy, Ilias Bilionis〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉State-of-the-art computer codes for simulating real physical systems are often characterized by vast number of input parameters. Performing uncertainty quantification (UQ) tasks with Monte Carlo (MC) methods is almost always infeasible because of the need to perform hundreds of thousands or even millions of forward model evaluations in order to obtain convergent statistics. One, thus, tries to construct a cheap-to-evaluate surrogate model to replace the forward model solver. For systems with large numbers of input parameters, one has to address the curse of dimensionality through suitable dimensionality reduction techniques. A popular class of dimensionality reduction methods are those that attempt to recover a low-dimensional representation of the high-dimensional feature space. However, such methods often tend to overestimate the intrinsic dimensionality of the input feature space. In this work, we demonstrate the use of deep neural networks (DNN) to construct surrogate models for numerical simulators. We parameterize the structure of the DNN in a manner that lends the DNN surrogate the interpretation of recovering a low-dimensional nonlinear manifold. The model response is a parameterized nonlinear function of the low-dimensional projections of the input. We think of this low-dimensional manifold as a nonlinear generalization of the notion of the 〈em〉active subspace〈/em〉. Our approach is demonstrated with a problem on uncertainty propagation in a stochastic elliptic partial differential equation (SPDE) with uncertain diffusion coefficient. We deviate from traditional formulations of the SPDE problem by lifting the assumption of fixed lengthscales of the uncertain diffusion field. Instead we attempt to solve a more challenging problem of learning a map between an arbitrary snapshot of the diffusion field and the response.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Jordi Casacuberta, Koen J. Groot, Henry J. Tol, Stefan Hickel〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Selective Frequency Damping (SFD) is a popular method for the computation of globally unstable steady-state solutions in fluid dynamics. The approach has two model parameters whose selection is generally unclear. In this article, a detailed analysis of the influence of these parameters is presented, answering several open questions with regard to the effectiveness, optimum efficiency and limitations of the method. In particular, we show that SFD is always capable of stabilising a globally unstable systems ruled by one unsteady unstable eigenmode and derive analytical formulas for optimum parameter values. We show that the numerical feasibility of the approach depends on the complex phase angle of the most unstable eigenvalue. A numerical technique for characterising the pertinent eigenmodes is presented. In combination with analytical expressions, this technique allows finding optimal parameters that minimise the spectral radius of a simulation, without having to perform an independent stability analysis. An extension to multiple unstable eigenmodes is derived. As computational example, a two-dimensional cylinder flow case is optimally stabilised using this method. We provide a physical interpretation of the stabilisation mechanism based on, but not limited to, this Navier–Stokes example.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): M. Esmaeilbeigi, O. Chatrabgoun, M. Cheraghi〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In the interpolation method, in some cases, one often has a number of data points and its derivatives, which are obtained by sampling or experimentation. In this case, the problem of finding an approximating function passing through these points and coinciding with given values of its derivatives at these points is generally known as “Hermite interpolation”. The Hermite interpolation is mostly a method of interpolating data points as a polynomial function that is faced with some challenges in high dimensions and on irregular domains. Radial basis functions take advantage of being flexible with respect to geometry, easy to implement in higher dimensions, and can also provide high order convergence. So, they can be applied as a suitable tool to high dimensional Hermite interpolation problem on irregular domains. In many applied systems, commonly available derivatives information is presented using fractional order derivatives instead of integer ones. For this purpose, in this paper, we assume that the values of an unknown function and its fractional derivatives at some distinct points are presented. Therefore, we intend to apply a new approach, which we call it as “fractional Hermite interpolation” with radial basis functions in high dimensions. Optimal recovery conditions for the fractional Hermite interpolant are investigated. Then, the existence and uniqueness of the solution in this type of generalized interpolation are proved. In order to increase the accuracy and stability of the method, Hilbert Schmidt's theory has also been used. The main advantages of the used method are its simplicity and efficiency in high dimensions, and over irregular domains. Finally, numerical results in one, two and three dimensions and a real-world problem are presented to support our theoretical analysis.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Majid Haghshenas, James A. Wilson, Ranganathan Kumar〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this work, we implement interfacial forces (surface tension, hydrostatic and viscous forces) by enlisting the finite volume discretization of GFM (Ghost Fluid Method) using A-CLSVOF (Algebraic Coupled Level-Set/Volume-of-Fluid) method for the mass conservative, and smooth interface description. The widely used PISO momentum solution to resolve the pressure–velocity coupling is presented along with the present GFM discretization and its placement within the PISO loop. The pressure jump at the interface due to the interfacial forces is made sharp via direct calculation of the modified pressure matrix coefficients corresponding to targeted interfacial cells, and as a source term for the jump value itself. The Level-Set field is enlisted for curvature computation in A-CLSVOF and for the interpolation and weighting of the relative contribution of the capillary force in adjacent for the matrix coefficients in the FV framework. To assess the A-CLSVOF/GFM performance, four canonical cases were studied. In the case of a static droplet in suspension, A-CLSVOF/GFM produces a sharp and accurate pressure jump compared to the traditional CSF implementation of A-CLSVOF. The interaction of viscous and capillary forces is proven to be accurate and consistent with theoretical results for the classical capillary wave. For the linear two-layer shear flow, GFM sharp treatment of the viscosity captured the velocity gradient across the interface and removed the diffusion of the viscous stresses caused by the discontinuous material properties. Finally, the combination of all GFM improvements proposed in this study are compared to experimental findings of terminal velocity for a gaseous bubble rising in a viscous fluid. GFM outperforms CSF with errors of 4.6% and 14.0% respectively.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Samira Nikkar, Jan Nordström〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉In this article, well-posedness and dual consistency of the linearized constant coefficient incompressible Navier–Stokes equations posed on time-dependent spatial domains are studied. To simplify the derivation of the dual problem and improve the accuracy of gradients, the second order formulation is transformed to first order form. Boundary conditions that simultaneously lead to boundedness of the primal and dual problems are derived.〈/p〉 〈p〉Fully discrete finite difference schemes on summation-by-parts form, in combination with the simultaneous approximation technique, are constructed. We prove energy stability and discrete dual consistency and show how to construct the penalty operators such that the scheme automatically adjusts to the variations of the spatial domain. As a result of the aforementioned formulations, stability and discrete dual consistency follow simultaneously.〈/p〉 〈p〉The method is illustrated by considering a deforming time-dependent spatial domain in two dimensions. The numerical calculations are performed using high order operators in space and time. The results corroborate the stability of the scheme and the accuracy of the solution. We also show that linear functionals are superconverging. Additionally, we investigate the convergence of non-linear functionals and the divergence of the solution.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Chris J. Budd, Andrew T.T. McRae, Colin J. Cotter〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In the context of numerical solution of PDEs, dynamic mesh redistribution methods (r-adaptive methods) are an important procedure for increasing the resolution in regions of interest, without modifying the connectivity of the mesh. Key to the success of these methods is that the mesh should be sufficiently refined (locally) and flexible in order to resolve evolving solution features, but at the same time not introduce errors through skewness and lack of regularity. Some state-of-the-art methods are bottom-up in that they attempt to prescribe both the local cell size and the alignment to features of the solution. However, the resulting problem is overdetermined, necessitating a compromise between these conflicting requirements. An alternative approach, described in this paper, is to prescribe only the local cell size and augment this an optimal transport condition to provide global regularity. This leads to a robust and flexible algorithm for generating meshes fitted to an evolving solution, with minimal need for tuning parameters. Of particular interest for geophysical modelling are meshes constructed on the surface of the sphere. The purpose of this paper is to demonstrate that meshes generated on the sphere using this optimal transport approach have good a-priori regularity and that the meshes produced are naturally aligned to various simple features. It is further shown that the sphere's intrinsic curvature leads to more regular meshes than the plane. In addition to these general results, we provide a wide range of examples relevant to practical applications, to showcase the behaviour of optimally transported meshes on the sphere. These range from axisymmetric cases that can be solved analytically to more general examples that are tackled numerically. Evaluation of the singular values and singular vectors of the mesh transformation provides a quantitative measure of the mesh anisotropy, and this is shown to match analytic predictions.〈/p〉〈/div〉 〈h5〉Graphical abstract〈/h5〉 〈div〉〈p〉〈/p〉〈/div〉 〈figure〉〈img src="https://ars.els-cdn.com/content/image/1-s2.0-S0021999118305515-gr001.jpg" width="494" alt="Graphical abstract for this article" title=""〉〈/figure〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Ahmad Al Takash, Marianne Beringhier, Mohammad Hammoud, Jean-Claude Grandidier〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Large computation time is widely considered to be the most important issue in scientific research especially in solving structural evolution problems. Recent developments in this domain have shown that the use of non-incremental schemes through Model Order Reduction led to important results in saving time. Yet, the question arises here how to obtain more time-saving. This paper examines an approach based on a collection of significant modes given by Proper Generalized Decomposition (PGD) solution for different time scales in order to save more computation time. The dictionary of the significant modes allows to construct an accurate solution for different characteristic times and different boundary problems compared to the full solution with a relative error rate less than 5% and with a large time saving of order 50 compared to Finite Element Method (FEM). The ability of the approach with respect to cycle time is discussed.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Christopher Lester, Christian A. Yates, Ruth E. Baker〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this work, we consider the problem of estimating summary statistics to characterise biochemical reaction networks of interest. Such networks are often described using the framework of the Chemical Master Equation (CME). For physically-realistic models, the CME is widely considered to be analytically intractable. A variety of Monte Carlo algorithms have therefore been developed to explore the dynamics of such networks empirically. Amongst them is the multi-level method, which uses estimates from multiple ensembles of sample paths of different accuracies to estimate a summary statistic of interest. In this work, we develop the multi-level method in two directions: (1) to increase the robustness, reliability and performance of the multi-level method, we implement an improved variance reduction method for generating the sample paths of each ensemble; and (2) to improve computational performance, we demonstrate the successful use of a different mechanism for choosing which ensembles should be included in the multi-level algorithm.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Dong-Yeop Na, Ben-Hur V. Borges, Fernando L. Teixeira〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present a finite-element time-domain (FETD) Maxwell solver for the analysis of body-of-revolution (BOR) geometries based on discrete exterior calculus (DEC) of differential forms and transformation optics (TO) concepts. We explore TO principles to map the original 3-D BOR problem to a 2-D one in the meridian 〈em〉ρz〈/em〉-plane based on a Cartesian coordinate system where the cylindrical metric is fully embedded into the constitutive properties of an effective inhomogeneous and anisotropic medium that fills the domain. The proposed solver uses a 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mo stretchy="false"〉(〈/mo〉〈msup〉〈mrow〉〈mtext〉TE〈/mtext〉〈/mrow〉〈mrow〉〈mi〉ϕ〈/mi〉〈/mrow〉〈/msup〉〈mo〉,〈/mo〉〈msup〉〈mrow〉〈mtext〉TM〈/mtext〉〈/mrow〉〈mrow〉〈mi〉ϕ〈/mi〉〈/mrow〉〈/msup〉〈mo stretchy="false"〉)〈/mo〉〈/math〉 field decomposition and an appropriate set of DEC-based basis functions on an irregular grid discretizing the meridian plane. A symplectic time discretization based on a leap-frog scheme is applied to obtain the full-discrete marching-on-time algorithm. We validate the algorithm by comparing the numerical results against analytical solutions for resonant fields in cylindrical cavities and against pseudo-analytical solutions for fields radiated by cylindrically symmetric antennas in layered media. We also illustrate the application of the algorithm for a particle-in-cell (PIC) simulation of beam-wave interactions inside a high-power backward-wave oscillator.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Vamsi Spandan, Detlef Lohse, Marco D. de Tullio, Roberto Verzicco〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this paper we propose and test the validity of simple and easy-to-implement algorithms within the immersed boundary framework geared towards large scale simulations involving thousands of deformable bodies in highly turbulent flows. First, we introduce a fast moving least squares (fast-MLS) approximation technique with which we speed up the process of building transfer functions during the simulations which leads to considerable reductions in computational time. We compare the accuracy of the fast-MLS against the exact moving least squares (MLS) for the standard problem of uniform flow over a sphere. In order to overcome the restrictions set by the resolution coupling of the Lagrangian and Eulerian meshes in this particular immersed boundary method, we present an adaptive Lagrangian mesh refinement procedure that is capable of drastically reducing the number of required nodes of the basic Lagrangian mesh when the immersed boundaries can move and deform. Finally, a coarse-grained collision detection algorithm is presented which can detect collision events between several Lagrangian markers residing on separate complex geometries with minimal computational overhead.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Peter Korn, Leonidas Linardakis〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A structure-preserving discretization of the shallow-water equations on unstructured spherical grids is introduced. The unstructured grids that we consider have triangular cells with a C-type staggering of variables, where scalar variables are located at centres of grid cells and normal components of velocity are placed at cell boundaries. The staggering necessitates reconstructions and these reconstructions are build into the algorithm such that the resulting discrete equations obey a weighted weak form. This approach, combined with a mimetic discretization of the differential operators of the shallow-water equations, provides a conservative discretization that preserves important aspects of the mathematical structure of the continuous equations, most notably the simultaneous conservation of quadratic invariants such as energy and enstrophy. The structure-preserving nature of our discretization is confirmed through theoretical analysis and through numerical experiments on two different triangular grids, a symmetrized icosahedral grid of nearly uniform resolution and a non-uniform triangular grid whose resolution increases towards the poles.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Louisa Schlachter, Florian Schneider〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Uncertainty Quantification through stochastic spectral methods is rising in popularity. We derive a modification of the classical stochastic Galerkin method, that ensures the hyperbolicity of the underlying hyperbolic system of partial differential equations. The modification is done using a suitable “slope” limiter, based on similar ideas in the context of kinetic moment models. We apply the resulting modified stochastic Galerkin method to the compressible Euler equations and the 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈msub〉〈mrow〉〈mtext〉M〈/mtext〉〈/mrow〉〈mrow〉〈mn〉1〈/mn〉〈/mrow〉〈/msub〉〈/math〉 model of radiative transfer. Our numerical results show that it can compete with other UQ methods like the intrusive polynomial moment method while being computationally inexpensive and easy to implement.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Navid Shervani-Tabar, Oleg V. Vasilyev〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This paper addresses one of the main challenges of the conservative level set method, namely the ill-conditioned behavior of the normal vector away from the interface. An alternative formulation for reconstruction of the interface is proposed. Unlike commonly used methods, which rely on a unit normal vector, the Stabilized Conservative Level Set (SCLS) makes use of a modified normal vector with diminishing magnitude away from the interface. With the new formulation, in the vicinity of the interface the reinitialization procedure utilizes compressive flux and diffusive terms only in normal direction with respect to the interface, thus, preserving the conservative level set properties, while away from the interface the directional diffusion mechanism automatically switches to homogeneous diffusion. The proposed formulation is robust and general. It is especially well suited for use with the adaptive mesh refinement (AMR) approaches, since for computational accuracy high resolution is only required in the vicinity of the interface, while away from the interface low resolution simulations might be sufficient. All of the results reported in this paper are obtained using the Adaptive Wavelet Collocation Method, a general arbitrary order AMR-type method, which utilizes wavelet decomposition to adapt on steep gradients in the solution while retaining a predetermined order of accuracy. Numerical solution for a number of benchmark problems has been carried out to demonstrate the performance of the SCLS method.〈/p〉〈/div〉