ALBERT

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〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): John W. Barrett, Harald Garcke, Robert Nürnberg〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Fourth order curvature driven interface evolution equations frequently appear in the natural sciences. Often axisymmetric geometries are of interest, and in this situation numerical computations are much more efficient. We will introduce and analyze several new finite element schemes for fourth order geometric evolution equations in an axisymmetric setting, and for selected schemes we will show existence, uniqueness and stability results. The presented schemes have very good mesh and stability properties, as will be demonstrated by several numerical examples.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Bryan D. Quaife, M. Nicholas J. Moore〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We develop numerical methods to simulate the fluid-mechanical erosion of many bodies in two-dimensional Stokes flow. The broad aim is to simulate the erosion of a porous medium (e.g. groundwater flow) with grain-scale resolution. Our fluid solver is based on a second-kind boundary integral formulation of the Stokes equations that is discretized with a spectrally-accurate Nyström method and solved with fast-multipole-accelerated GMRES. The fluid solver provides the surface shear stress which is used to advance solid boundaries. We regularize interface evolution via curvature penalization using the 〈em〉θ〈/em〉–〈em〉L〈/em〉 formulation, which affords numerically stable treatment of stiff terms and therefore permits large time steps. The overall accuracy of our method is spectral in space and second-order in time. The method is computationally efficient, with the fluid solver requiring 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mi mathvariant="script"〉O〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈mi〉N〈/mi〉〈mo stretchy="false"〉)〈/mo〉〈/math〉 operations per GMRES iteration, a mesh-independent number of GMRES iterations, and a one-time 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mi mathvariant="script"〉O〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈msup〉〈mrow〉〈mi〉N〈/mi〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈/mrow〉〈/msup〉〈mo stretchy="false"〉)〈/mo〉〈/math〉 computation to compute the shear stress. We benchmark single-body results against analytical predictions for the limiting morphology and vanishing rate. Multibody simulations reveal the spontaneous formation of channels between bodies of close initial proximity. The channelization is associated with a dramatic reduction in the resistance of the porous medium, much more than would be expected from the reduction in grain size alone.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 378〈/p〉 〈p〉Author(s): L. Grosheintz-Laval, R. Käppeli〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A high-order well-balanced scheme for the Euler equations with gravitation is presented. The scheme is able to preserve a spatially high-order accurate discrete representation of isentropic hydrostatic equilibria. It is based on a novel local hydrostatic reconstruction, which, in combination with any standard high-order accurate reconstruction procedure, achieves genuine high-order accuracy for smooth solutions close or away from equilibrium. The resulting scheme is very simple and can be implemented into any existing finite volume code with minimal effort. Moreover, the scheme is not tied to any particular form of the equation of state, which is crucial for example in astrophysical applications. Several numerical experiments were performed with a third-order accurate reconstruction. They demonstrate the robustness and high-order accuracy of the scheme nearby and out of hydrostatic equilibrium.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 378〈/p〉 〈p〉Author(s): Shashank Jaiswal, Alina A. Alexeenko, Jingwei Hu〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The Boltzmann equation, an integro-differential equation for the molecular distribution function in the physical and velocity phase space, governs the fluid flow behavior at a wide range of physical conditions, including compressible, turbulent, as well as flows involving further physics such as non-equilibrium internal energy exchange and chemical reactions. Despite its wide applicability, deterministic solution of the Boltzmann equation presents a huge computational challenge, and often the collision operator is simplified for practical reasons. In this work, we introduce a highly accurate deterministic method for the full Boltzmann equation which couples the Runge–Kutta discontinuous Galerkin (RKDG) discretization in time and physical space (Su et al. (2015) [1]) and the recently developed fast Fourier spectral method in velocity space (Gamba et al. (2017) [2]). The novelty of this approach encompasses three aspects: first, the fast spectral method for the collision operator applies to general collision kernels with little or no practical limitations, and in order to adapt to the spatial discretization, we propose here a singular-value-decomposition based algorithm to further reduce the cost in evaluating the collision term; second, the DG formulation employed has high order of accuracy at element-level, and has shown to be more efficient than the finite volume method; thirdly, the element-local compact nature of DG as well as our collision algorithm is amenable to effective parallelization on massively parallel architectures. The solver has been verified against analytical Bobylev–Krook–Wu solution. Further, the standard benchmark test cases of rarefied Fourier heat transfer, Couette flow, oscillatory Couette flow, normal shock wave, lid-driven cavity flow, and thermally driven cavity flow have been studied and their results are compared against direct simulation Monte Carlo (DSMC) solutions with equivalent molecular collision models or published deterministic solutions.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 378〈/p〉 〈p〉Author(s): Kossi-Mensah Kodjo, Julien Yvonnet, Mustapha Karkri, Karam Sab〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this work, a multiscale model for thermomechanical properties of composite structures containing phase change particles is developed. For the mechanical part, a classical linear computational homogenization procedure is employed. For the thermal part, due to the strong nonlinear, history-dependent thermal effects, a concurrent multiscale (FE〈sup〉2〈/sup〉) method is extended to take into account the presence of Phase Change Materials particles (PCM) at the microscale. The PCM inclusions change from liquid to solid state in the range of room temperature. This phase change induces a modified macroscopic thermal behavior, which can be used e.g. to design materials with enhanced thermal inertia and reduce energy consumption in civil engineering constructions. The technique allows taking into account accurately the fully nonlinear, history-dependent thermal behavior through numerical calculations at the microscale based on a Representative Volume Element (RVE) and its effect at the macroscale. The method is applied to concrete material including paraffin wax PCM. The results show the benefits of the PCM on the thermal behavior, including shifted and smoothed temperature response as compared to materials without PCM particles.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 378〈/p〉 〈p〉Author(s): Andrés M. Rueda-Ramírez, Juan Manzanero, Esteban Ferrer, Gonzalo Rubio, Eusebio Valero〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉High-order discontinuous Galerkin methods have become a popular technique in computational fluid dynamics because their accuracy increases spectrally in smooth solutions with the order of the approximation. However, their main drawback is that increasing the order also increases the computational cost. Several techniques have been introduced in the past to reduce this cost. On the one hand, local mesh adaptation strategies based on error estimation have been proposed to reduce the number of degrees of freedom while keeping a similar accuracy. On the other hand, multigrid solvers may accelerate time marching computations for a fixed number of degrees of freedom.〈/p〉 〈p〉In this paper, we combine both methods and present a novel anisotropic p-adaptation multigrid algorithm for steady-state problems that uses the multigrid scheme both as a solver and as an anisotropic error estimator. To achieve this, we show that a recently developed anisotropic truncation error estimator [1, A. M. Rueda-Ramírez, G. Rubio, E. Ferrer, E. Valero, Truncation error estimation in the p-anisotropic discontinuous Galerkin spectral element method, J. Sci. Comput.] is perfectly suited to be performed inside the multigrid cycle with negligible extra cost. Furthermore, we introduce a multi-stage p-adaptation procedure which reduces the computational time when very accurate results are required.〈/p〉 〈p〉The proposed methods are tested for the compressible Navier–Stokes equations, where we investigate two steady-state cases. First, the 2D boundary layer flow on a flat plate is studied to assess accuracy and computational cost of the algorithm, where a speed-up of 816 is achieved compared to the traditional explicit method. Second, the 3D flow around a sphere is simulated and used to test the anisotropic properties of the proposed method, where a speed-up of 152 is achieved compared to the explicit method. The proposed multi-stage procedure achieved a speed-up of 2.6 in comparison to the single-stage method in highly accurate simulations.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): Ryan Galagusz, Steve McFee〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The purpose of this research is to describe an efficient iterative method suitable for obtaining high accuracy solutions to high frequency time-harmonic scattering problems. The method allows for both refinement of local polynomial degree and non-conforming mesh refinement, including multiple hanging nodes per edge. Rather than globally assemble the finite element system, we describe an iterative domain decomposition method which can use element-wise fast solvers for elements of large degree. Since continuity between elements is enforced through moment equations, the resulting constraint equations are hierarchical. We show that, for high frequency problems, a subset of these constraints should be directly enforced, providing the coarse space in the dual-primal domain decomposition method. The subset of constraints is chosen based on a dispersion criterion involving mesh size and wavenumber. By increasing the size of the coarse space according to this criterion, the number of iterations in the domain decomposition method depends only weakly on the wavenumber. We demonstrate this convergence behaviour on standard domain decomposition test problems and conclude the paper with application of the method to electromagnetic problems in two dimensions. These examples include beam steering by lenses and photonic crystal waveguides, as well as radar cross section computation for dielectric, perfect electric conductor, and electromagnetic cloak scatterers.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 378〈/p〉 〈p〉Author(s): Jesse Chan, Lucas C. Wilcox〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We construct entropy conservative and entropy stable discontinuous Galerkin (DG) discretizations for time-dependent nonlinear hyperbolic conservation laws on curvilinear meshes. The resulting schemes preserve a semi-discrete quadrature approximation of a continuous global entropy inequality. The proof requires the satisfaction of a discrete geometric conservation law, which we enforce through an appropriate polynomial approximation. We extend the construction of entropy conservative and entropy stable DG schemes to the case when high order curvilinear mass matrices are approximated using low-storage weight-adjusted approximations, and describe how to retain global conservation properties under such an approximation. For certain types of curvilinear meshes, these weight-adjusted approximations deliver optimal rates of convergence. Finally, the high order accuracy, local conservation, and discrete conservation or dissipation of entropy for under-resolved solutions are verified through numerical experiments for the compressible Euler equations on triangular and tetrahedral meshes.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 378〈/p〉 〈p〉Author(s): M. Jedouaa, C.-H. Bruneau, E. Maitre〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉An efficient method to capture an arbitrary number of fluid/structure interfaces in a level-set framework is built, following ideas introduced for contour capturing in image analysis. Using only three label maps and two distance functions we succeed in locating and evolving the bodies independently in the whole domain and get the distance between the closest bodies in order to apply a collision force whatever the number of cells is. The method is applied to rigid solid bodies in order to compare to the results available in the literature. In that case, a global penalization model uses the label maps to follow the solid bodies all together without a separate computation of each body velocity. Numerical simulations are performed in two- and three-dimensions. An application to immersed vesicles is also proposed and shows the capability and efficiency of the method to handle numerical contacts between elastic bodies at low resolution. Two-dimensional simulations of vesicles under various flow conditions are presented.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 378〈/p〉 〈p〉Author(s): Wei Su, Peng Wang, Haihu Liu, Lei Wu〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The Couette flow is one of the fundamental problems of rarefied gas dynamics, which has been investigated extensively based on the linearized Boltzmann equation (LBE) of hard-sphere molecules and simplified kinetic model equations. However, how the different intermolecular potentials affect the viscous slip coefficient and the structure of Knudsen layer remains unclear. Here, a novel synthetic iteration scheme (SIS) is developed for the LBE to find solutions of Couette flow accurately and efficiently: the velocity distribution function is first solved by the conventional iterative scheme, then it is modified such that in each iteration i) the flow velocity is guided by an ordinary differential equation that is asymptotic-preserving at the Navier–Stokes limit and ii) the shear stress is equal to the average shear stress. Based on the Bhatnagar–Gross–Krook model, the SIS is assessed to be efficient and accurate. Then we investigate the Knudsen layer function for gases interacting through the inverse power-law, shielded Coulomb, and Lennard-Jones potentials, subject to diffuse-specular and Cercignani–Lampis gas-surface boundary conditions. When the tangential momentum accommodation coefficient (TMAC) is not larger than one, the Knudsen layer function is strongly affected by the potential, where its value and width increase with the effective viscosity index of gas molecules. Moreover, the Knudsen layer function exhibits similarities among different values of TMAC when the intermolecular potential is fixed. For Cercignani–Lampis boundary condition with TMAC larger than one, both the viscous slip coefficient and Knudsen layer function are affected by the intermolecular potential, especially when the “backward” scattering limit is approached. With the asymptotic theory by Jiang and Luo (2016) [14] for the singular behavior of the velocity gradient in the vicinity of solid surfaces, we find that the whole Knudsen layer function can be well fitted by the power series 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈msubsup〉〈mrow〉〈mo〉∑〈/mo〉〈/mrow〉〈mrow〉〈mi〉n〈/mi〉〈mo〉=〈/mo〉〈mn〉0〈/mn〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈/mrow〉〈/msubsup〉〈msubsup〉〈mrow〉〈mo〉∑〈/mo〉〈/mrow〉〈mrow〉〈mi〉m〈/mi〉〈mo〉=〈/mo〉〈mn〉0〈/mn〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈/mrow〉〈/msubsup〉〈msub〉〈mrow〉〈mi〉c〈/mi〉〈/mrow〉〈mrow〉〈mi〉n〈/mi〉〈mo〉,〈/mo〉〈mi〉m〈/mi〉〈/mrow〉〈/msub〉〈msup〉〈mrow〉〈mi〉x〈/mi〉〈/mrow〉〈mrow〉〈mi〉n〈/mi〉〈/mrow〉〈/msup〉〈msup〉〈mrow〉〈mo stretchy="false"〉(〈/mo〉〈mi〉x〈/mi〉〈mi mathvariant="normal"〉ln〈/mi〉〈mo〉⁡〈/mo〉〈mi〉x〈/mi〉〈mo stretchy="false"〉)〈/mo〉〈/mrow〉〈mrow〉〈mi〉m〈/mi〉〈/mrow〉〈/msup〉〈/math〉, where 〈em〉x〈/em〉 is the distance to the solid surface. Finally, the experimental data of the Knudsen layer profile are explained by the LBE solution with proper values of the viscosity index and TMAC.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 378〈/p〉 〈p〉Author(s): Alessandro Franci, Massimiliano Cremonesi〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This paper proposes two regularized models of the 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mi〉μ〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈mi〉I〈/mi〉〈mo stretchy="false"〉)〈/mo〉〈/math〉-rheology and shows their application to the numerical simulation of 3D dense granular flows. The proposed regularizations are inspired by the Papanastasiou and Bercovier–Engleman methods, typically used to approximate the Bingham law. The key idea is to keep limited the value of the apparent viscosity for low shear rates without introducing a fixed cutoff. The proposed techniques are introduced into the Particle Finite Element Method (PFEM) framework to deal with the large deformations expected in free-surface granular flows. After showing the numerical drawbacks associated to the standard 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mi〉μ〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈mi〉I〈/mi〉〈mo stretchy="false"〉)〈/mo〉〈/math〉-rheology, the two regularization strategies are derived and discussed. The regularized 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mi〉μ〈/mi〉〈mo stretchy="false"〉(〈/mo〉〈mi〉I〈/mi〉〈mo stretchy="false"〉)〈/mo〉〈/math〉-rheology is then applied to the simulation of the collapse of 2D and 3D granular columns. The numerical results show that the regularization techniques improve substantially the conditioning of the linear system without affecting the solution accuracy. A good agreement with the experimental tests and other numerical methods is obtained in all the analyzed problems.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): Thomas Adams, Stefano Giani, William M. Coombs〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this paper, an efficient, high-order accurate, level set reinitialisation method is proposed, based on the elliptic reinitialisation method (Basting and Kuzmin, 2013 [1]), which is discretised spatially using the discontinuous Galerkin (DG) symmetric interior penalty method (SIPG). In order to achieve this a number of improvements have been made to the elliptic reinitialisation method including; reformulation of the underlying minimisation problem driving the solution; adoption of a Lagrange multiplier approach for enforcing a Dirichlet boundary condition on the implicit level set interface; and adoption of a narrow band approach. Numerical examples confirm the high-order accuracy of the resultant method by demonstrating experimental orders of convergence congruent with optimal convergence rates for the SIPG method, that is 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈msup〉〈mrow〉〈mi〉h〈/mi〉〈/mrow〉〈mrow〉〈mi〉p〈/mi〉〈mo〉+〈/mo〉〈mn〉1〈/mn〉〈/mrow〉〈/msup〉〈/math〉 and 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈msup〉〈mrow〉〈mi〉h〈/mi〉〈/mrow〉〈mrow〉〈mi〉p〈/mi〉〈/mrow〉〈/msup〉〈/math〉 in the 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si3.gif" overflow="scroll"〉〈msup〉〈mrow〉〈mi〉L〈/mi〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈/mrow〉〈/msup〉〈/math〉 and DG norms respectively. Furthermore, the degree to which the level set function satisfies the Eikonal equation improves proportionally to 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈msup〉〈mrow〉〈mi〉h〈/mi〉〈/mrow〉〈mrow〉〈mi〉p〈/mi〉〈/mrow〉〈/msup〉〈/math〉, and the often ignored homogeneous Dirichlet boundary condition on the interface is shown to be satisfied accurately with a rate of convergence of at least 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si4.gif" overflow="scroll"〉〈msup〉〈mrow〉〈mi〉h〈/mi〉〈/mrow〉〈mrow〉〈mn〉2〈/mn〉〈/mrow〉〈/msup〉〈/math〉 for all polynomial orders.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 378〈/p〉 〈p〉Author(s): 〈/p〉

Description: 〈p〉Publication date: Available online 19 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): Benjamin Aymard, Urbain Vaes, Marc Pradas, Serafim Kalliadasis〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We propose a new numerical method to solve the Cahn-Hilliard equation coupled with non-linear wetting boundary conditions. We show that the method is mass-conservative and that the discrete solution satisfies a discrete energy law similar to the one satisfied by the exact solution. We perform several tests inspired by realistic situations to verify the accuracy and performance of the method: wetting of a chemically heterogeneous substrate in three dimensions, wetting-driven nucleation in a complex two-dimensional domain and three-dimensional diffusion through a porous medium.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): Zachary J. Silberman, Thomas R. Adams, Joshua A. Faber, Zachariah B. Etienne, Ian Ruchlin〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Many codes have been developed to study highly relativistic, magnetized flows around and inside compact objects. Depending on the adopted formalism, some of these codes evolve the vector potential 〈strong〉A〈/strong〉, and others evolve the magnetic field 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mi mathvariant="bold"〉B〈/mi〉〈mo〉=〈/mo〉〈mi mathvariant="normal"〉∇〈/mi〉〈mo〉×〈/mo〉〈mi mathvariant="bold"〉A〈/mi〉〈/math〉 directly. Given that these codes possess unique strengths, sometimes it is desirable to start a simulation using a code that evolves 〈strong〉B〈/strong〉 and complete it using a code that evolves 〈strong〉A〈/strong〉. Thus transferring the data from one code to another would require an inverse curl algorithm. This paper describes two new inverse curl techniques in the context of Cartesian numerical grids: a cell-by-cell method, which scales approximately linearly with the numerical grid, and a global linear algebra approach, which has worse scaling properties but is generally more robust (e.g., in the context of a magnetic field possessing some nonzero divergence). We demonstrate these algorithms successfully generate smooth vector potential configurations in challenging special and general relativistic contexts.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): Shuo Yang, Samuel F. Potter, Maria K. Cameron〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Nongradient SDEs with small white noise often arise when modeling biological and ecological time-irreversible processes. If the governing SDE were gradient, the maximum likelihood transition paths, transition rates, expected exit times, and the invariant probability distribution would be given in terms of its potential function. The quasipotential plays a similar role for nongradient SDEs. Unfortunately, the quasipotential is the solution of a functional minimization problem that can be obtained analytically only in some special cases. We propose a Dijkstra-like solver for computing the quasipotential on regular rectangular meshes in 3D. This solver results from a promotion and an upgrade of the previously introduced ordered line integral method with the midpoint quadrature rule for 2D SDEs. The key innovations that have allowed us to keep the CPU times reasonable while maintaining good accuracy are (〈em〉i〈/em〉) a new hierarchical update strategy, 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mo stretchy="false"〉(〈/mo〉〈mi〉i〈/mi〉〈mi〉i〈/mi〉〈mo stretchy="false"〉)〈/mo〉〈/math〉 the use of Karush–Kuhn–Tucker theory for rejecting unnecessary simplex updates, and 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mo stretchy="false"〉(〈/mo〉〈mi〉i〈/mi〉〈mi〉i〈/mi〉〈mi〉i〈/mi〉〈mo stretchy="false"〉)〈/mo〉〈/math〉 pruning the number of admissible simplexes and a fast search for them. An extensive numerical study is conducted on a series of linear and nonlinear examples where the quasipotential is analytically available or can be found at transition states by other methods. In particular, the proposed solver is applied to Tao's examples where the transition states are hyperbolic periodic orbits, and to a genetic switch model by Lv et al. (2014) [21]. The C source code implementing the proposed algorithm is available at M. Cameron's web page.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): Baoli Yin, Yang Liu, Hong Li, Siriguleng He〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this article, a fast algorithm based on time two-mesh (TT-M) finite element (FE) scheme, which aims at solving nonlinear problems quickly, is considered to numerically solve the nonlinear space fractional Allen–Cahn equations with smooth and non-smooth solutions. The implicit second-order 〈em〉θ〈/em〉 scheme containing both implicit Crank–Nicolson scheme and second-order backward difference method is applied to time direction, a fast TT-M method is used to increase the speed of calculation, and the FE method is developed to approximate the spacial direction. The TT-M FE algorithm includes the following main computing steps: firstly, a nonlinear implicit second-order 〈em〉θ〈/em〉 FE scheme on the time coarse mesh 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈msub〉〈mrow〉〈mi〉τ〈/mi〉〈/mrow〉〈mrow〉〈mi〉c〈/mi〉〈/mrow〉〈/msub〉〈/math〉 is solved by a nonlinear iterative method; secondly, based on the chosen initial iterative value, a linearized FE system on time fine mesh 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mi〉τ〈/mi〉〈mo〉〈〈/mo〉〈msub〉〈mrow〉〈mi〉τ〈/mi〉〈/mrow〉〈mrow〉〈mi〉c〈/mi〉〈/mrow〉〈/msub〉〈/math〉 is solved, where some useful coarse numerical solutions are found by the Lagrange's interpolation formula. The analysis for both stability and a priori error estimates are made in detail. Finally, three numerical examples with smooth and non-smooth solutions are provided to illustrate the computational efficiency in solving nonlinear partial differential equations, from which it is easy to find that the computing time can be saved.〈/p〉〈/div〉

Description: 〈p〉Publication date: Available online 13 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics〈/p〉 〈p〉Author(s): X.-X. Cai, T. Kittelmann, E. Klinkby, J.I. Márquez Damián〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Distributions of inelastically scattered neutrons can be quantum dynamically described by a scattering kernel. We present an accurate and computationally efficient rejection method for sampling a given scattering kernel of any isotropic material. The proposed method produces continuous neutron energy and angular distributions, typically using just a single interpolation per sampling. We benchmark the results of this method against those from accurate analytical models and one of the major neutron transport codes. We also show the results of applying this method to the conventional discrete double differential cross sections.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): Hailong Guo, Xu Yang, Yi Zhu〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Photonic graphene, a photonic crystal with honeycomb structures, has been intensively studied in both theoretical and applied fields. Similar to graphene which admits Dirac Fermions and topological edge states, photonic graphene supports novel and subtle propagating modes (edge modes) of electromagnetic waves. These modes have wide applications in many optical systems. In this paper, we propose a novel gradient recovery method based on Bloch theory for the computation of topological edge modes in photonic graphene. Compared to standard finite element methods, this method provides higher order accuracy with the help of gradient recovery technique. This high order accuracy is desired for constructing the propagating electromagnetic modes in applications. We analyze the accuracy and prove the superconvergence of this method. Numerical examples are presented to show the efficiency by computing the edge mode for the 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mi mathvariant="script"〉P〈/mi〉〈/math〉-symmetry and 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mi mathvariant="script"〉C〈/mi〉〈/math〉-symmetry breaking cases in honeycomb structures.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): Andrew Christlieb, Wei Guo, Yan Jiang〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this paper, a class of high order numerical schemes is proposed for solving Hamilton–Jacobi (H–J) equations. This work is regarded as an extension of our previous work for nonlinear degenerate parabolic equations, see Christlieb et al. [14], which relies on a special kernel-based formulation of the solutions and successive convolution. When applied to the H–J equations, the newly proposed scheme attains genuinely high order accuracy in both space and time, and more importantly, it is unconditionally stable, hence allowing for much larger time step evolution compared with other explicit schemes and saving computational cost. A high order weighted essentially non-oscillatory methodology and a novel nonlinear filter are further incorporated to capture the correct viscosity solution. Furthermore, by coupling the recently proposed inverse Lax–Wendroff boundary treatment technique, this method is very flexible in handing complex geometry as well as general boundary conditions. We perform numerical experiments on a collection of numerical examples, including H–J equations with linear, nonlinear, convex or non-convex Hamiltonians. The efficacy and efficiency of the proposed scheme in approximating the viscosity solution of general H–J equations is verified.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 February 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 379〈/p〉 〈p〉Author(s): David Luquet, Régis Marchiano, François Coulouvrat〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Acoustical shock waves can be generated by numerous atmospheric sources, either natural – like thunder and volcanoes – or anthropic – like explosions, sonic boom or buzz saw noise. The prediction of their long-range propagation remains a numerical challenge at 3D because of the large propagation distance to wavelength ratio, and of the high frequency / small wavelength content associated to shocks. In this paper, an original numerical method for propagating acoustical shock waves in three-dimensional heterogeneous media is proposed. Heterogeneities can result from temperature or density gradients and also from atmospheric shear and turbulent flows. The method called FLHOWARD (for FLow and Heterogeneities in a One-Way Approximation of the nonlineaR wave equation in 3D) is based on a one-way solution of a generalized nonlinear wave equation. Even though backscattering is neglected, it does not suffer from the limitations of classical ray theory nor from the angular limitations of the popular parabolic methods. The numerical approach is based on a split-step method, which has the advantage of splitting the original equation into simpler ones associated with specific physical mechanisms: diffraction, flows, heterogeneities, nonlinearities, absorption and relaxation. The method has been developed on parallel architecture for very high demanding 3D configurations using the Single Method Multiple Data paradigm. The method is validated through several test cases. A study of the lateral cut-off of the sonic boom finally illustrates the potentialities of the method for realistic cases.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Darren Engwirda〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The generation of high-quality staggered unstructured grids is considered, leading to the development of a new optimisation-based strategy designed to construct weighted ‘Regular-Power’ tessellations appropriate for co-volume type numerical discretisation techniques. This new framework aims to extend the conventional Delaunay–Voronoi primal-dual structure; seeking to assemble generalised orthogonal tessellations with enhanced geometric quality. The construction of these grids is motivated by the desire to improve the performance and accuracy of numerical methods based on unstructured co-volume type schemes, including various staggered grid techniques for the simulation of fluid dynamics and hyperbolic transport. In this study, a new hybrid optimisation strategy is proposed; seeking to optimise the geometry, topology and weights associated with general, two-dimensional Regular-Power tessellations using a combination of gradient-ascent and energy-based techniques. The performance of this new method is tested experimentally, with a range of complex, multi-resolution primal-dual grids generated for various coastal and regional ocean modelling applications.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): R. Aubry, B.K. Karamete, E.L. Mestreau, C. Jones, S. Dey〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉Two particular aspects of volume boundary layer mesh generation applied to non smooth geometries are considered in this work. First, the treatment of concave ridges and corners is tackled from a generic viewpoint. Entropy satisfying elements are generated where shocks form in the volume. This proves useful to avoid premature halt of the boundary layer, and therefore potential jumps in the normal size. The connection with the Voronoi diagram is commented. Second, boundary layer adaptivity in the tangential plane is considered to honor arbitrary sizing prescription, and avoid size mismatch between the boundary layer and the isotropic sizing.〈/p〉 〈p〉It is shown that a strict semi-structured framework has to be abandoned in general to accommodate changes in the mesh topology. Size transition between boundary layer and fully unstructured anisotropic mesh is automatically taken into account. Both the concavity problem and the tangential adaptivity are presented together, since they require similar mesh operators. Various numerical examples illustrate the method.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Wen Yan, Michael Shelley〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉It is well-known that by placing judiciously chosen image point forces and doublets to the Stokeslet above a flat wall, the no-slip boundary condition can be conveniently imposed on the wall Blake (1971) [8]. However, to further impose periodic boundary conditions on directions parallel to the wall usually involves tedious derivations because single or double periodicity in Stokes flow may require the periodic unit to have no net force, which is not satisfied by the well-known image system. In this work we present a force-neutral image system. This neutrality allows us to represent the Stokes image system in a universal formulation for non-periodic, singly periodic and doubly periodic geometries. This formulation enables the black-box style usage of fast kernel summation methods. We demonstrate the efficiency and accuracy of this new image method with the periodic kernel independent fast multipole method in both non-periodic and periodic geometries. We then extend this new image system to other widely used Stokes fundamental solutions, including the Laplacian of the Stokeslet and the Rotne–Prager–Yamakawa tensor.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Natalie N. Beams, Andreas Klöckner, Luke N. Olson〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This paper presents a high-order method for solving an interface problem for the Poisson equation on embedded meshes through a coupled finite element and integral equation approach. The method is capable of handling homogeneous or inhomogeneous jump conditions without modification and retains high-order convergence close to the embedded interface. We present finite element–integral equation (FE–IE) formulations for interior, exterior, and interface problems. The treatments of the exterior and interface problems are new. The resulting linear systems are solved through an iterative approach exploiting the second-kind nature of the IE operator combined with algebraic multigrid preconditioning for the FE part. Assuming smooth continuations of coefficients and right-hand-side data, we show error analysis supporting high-order accuracy. Numerical evidence further supports our claims of efficiency and high-order accuracy for smooth data.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Philippe Chartier, Nicolas Crouseilles, Xiaofei Zhao〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this paper, we consider numerical methods for solving the two-dimensional Vlasov–Poisson equation in the finite Larmor radius approximation regime. The model describes the behaviour of charged particles under a strong external magnetic field and the finite Larmor radius scaling. We discretise the equation under Particle-in-Cell method, where the characteristics equations are highly oscillatory system in the limit regime. We apply popular numerical integrators including splitting methods, multi-revolution composition methods, two-scale formulation method and limit solver to integrate the characteristics. We then highlight the strengths and drawbacks of each method. Finally, numerical experiments are presented, and comparisons on the accuracy, efficiency and long-time behaviour of the methods are made, pointing to the method with the best performance for this problem.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Hojun You, Chongam Kim〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The present paper deals with a new improvement of hierarchical multi-dimensional limiting process for resolving the subcell distribution of high-order methods on two-dimensional mixed meshes. From previous studies, the multi-dimensional limiting process (MLP) was hierarchically extended to the discontinuous Galerkin (DG) method and the flux reconstruction/correction procedure via reconstruction (FR/CPR) method on simplex meshes. It was reported that the hierarchical MLP (〈em〉h〈/em〉MLP) shows several remarkable characteristics such as the preservation of the formal order-of-accuracy in smooth region and a sharp capturing of discontinuities in an efficient and accurate manner. At the same time, it was also surfaced that such characteristics are valid only on simplex meshes, and numerical Gibbs–Wilbraham oscillations are concealed in subcell distribution in the form of high-order polynomial modes. Subcell Gibbs–Wilbraham oscillations become potentially unstable near discontinuities and adversely affect numerical solutions in the sense of cell-averaged solutions as well as subcell distributions. In order to overcome the two issues, the behavior of the 〈em〉h〈/em〉MLP on mixed meshes is mathematically examined, and the simplex-decomposed 〈em〉P1〈/em〉-projected MLP condition and smooth extrema detector are derived. Secondly, a troubled-boundary detector is designed by analyzing the behavior of computed solutions across boundary-edges. Finally, 〈em〉h〈/em〉MLP_BD is proposed by combining the simplex-decomposed 〈em〉P1〈/em〉-projected MLP condition and smooth extrema detector with the troubled-boundary detector. Through extensive numerical tests, it is confirmed that the 〈em〉h〈/em〉MLP_BD scheme successfully eliminates subcell oscillations and provides reliable subcell distributions on two-dimensional triangular grids as well as mixed grids, while preserving the expected order-of-accuracy in smooth region.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Quentin Carmouze, François Fraysse, Richard Saurel, Boniface Nkonga〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A simple method is developed to couple accurately translational motion of rigid bodies to compressible fluid flows. Solid rigid bodies are tracked through a Level-Set function. Numerical diffusion is controlled thanks to a compressive limiter (Overbee) in the frame of MUSCL-type-scheme, giving an excellent compromise between accuracy and efficiency on unstructured meshes [9]. The method requires low resolution to preserve solid bodies' volume. Several coupling methods are then addressed to couple rigid body motion to fluid flow dynamics: a method based on stiff relaxation and two methods based on Ghost cells [13] and immersed boundaries. Their accuracy and convergence rates are compared against an immersed piston problem in 1D having exact solution. The second Ghost cell method is shown to be the most efficient. It is then extended to multidimensional computations on unstructured meshes and its accuracy is checked against flow computations around cylindrical bodies. Reference results are obtained when the flow evolves around a rigid body at rest. The same rigid body is then considered with prescribed velocity moving in a flow at rest. Computed results involving wave dynamics match very well. The method is then extended to two-way coupling and illustrated to several examples involving shock wave interaction with solid particles as well as phase transition induced by projectiles motion in liquid–gas mixtures.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Jie Zhang, Ming-Jiu Ni〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A new phase change model has been developed for the simulation of incompressible multiphase magnetohydrodynamics based on the Volume-of-Fluid method. To decrease the pressure oscillations when large density contrasts are present between the liquid phase and the vapor phase, a smooth distribution of sharp mass transfer rate within a narrow region surrounding the interface is adopted, and a ghost-cell approach is used to impose the saturating temperature at the liquid–vapor interface when solving the energy equation. After that, the method has been implemented in an incompressible multiphase magnetohydrodynamics solver developed in our previous work (Zhang and Ni (2014) [3]). Moreover, when computing the electromagnetic fields, a cut-cell approach is implemented to keep the sharpness of the interface, which is treated as an electrically insulating boundary as it translates and deforms with the fluid. The phase change model has been verified for a series of one-dimensional, two-dimensional and three-dimensional problems, while the numerical results agree well with either the theoretical solutions or the experimental data. In particular, by simulating the vapor bubble rising in superheated liquid under nonzero gravity in presence of external magnetic field, the magnetohydrodynamics effect on the vaporization of the rising bubble is investigated and we observe the magnetic fields to suppress the vapor bubble growth during the phase change. At last, both two-dimensional and three-dimensional film boiling simulations are conducted, which show good qualitative agreement with heat transfer correlations, and the vapor bubble is found to elongate along the direction of the magnetic field during its growth, moreover, the time instant for the vapor bubble to detach from the film is also delayed.〈/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. Frantzis, D.G.E. Grigoriadis〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉The numerical simulation of two-fluid flows with sharp interfaces is a challenging field, not only because of their complicated physical mechanisms, but also because of increased computational cost. An efficient and robust numerical formulation for incompressible two-fluid flows is proposed. Its novelty is the consistent coupling of Fast Direct Solvers (FDS) with the Immersed Boundary (IB) method to represent solid boundaries. Such a coupling offers several advantages. First, it extends the range of applicability of the IB method. Second, it allows the simulation of practical problems in geometrically complicated domains at a significantly reduced cost. Third, it can shed light on regions of the parametric space which are considered out of reach, or even impossible today.〈/p〉 〈p〉Instead of using a conventional variable coefficient pressure Poisson equation, a pressure-correction scheme is suggested for the solution of a constant coefficient Poisson equation for the pressure difference, extending the novel work of Dodd and Ferrante [8]. The conservative Level-set (LS) method is used to track the interface between the two fluids. Appropriate schemes, based on the local directional Ghost Cell Approach (GCA) are proposed, in order to satisfy the boundary conditions (BCs) of the pressure and the LS function around the IB.〈/p〉 〈p〉The accuracy, robustness, and performance of the proposed method is demonstrated by several validations against conventional approaches and experiments. The results verify that the pressure BCs are properly recovered along the IB solid interface, while a non-smooth pressure field is also allowed across the solid obstacle. The accuracy of the method was found to be 2nd-order, both in time and space. The performance of the proposed method is compared against the conventional approach using a multigrid iterative solver. The impact of the time-step on the accuracy of the constant coefficient approach is examined. Results show that the final speed-up strongly depends on the specific physical and numerical parameters such as the density ratio or the Reynolds number. It is demonstrated that for the range of parameters examined, speed-up factors of 100–10 can be achieved for density ratios of 10–1000 respectively.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): W. Tierens〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this paper we introduce an explicit and provably conditionally stable Finite Difference Time Domain (FDTD) algorithm for Maxwell's equations, with local refinement in both the spatial discretization length and in the time step (spatiotemporal refinement). This enables local spatial refinement with a locally reduced time step.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Kun Luo, Zhuo Wang, Junhua Tan, Jianren Fan〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The direct-forcing immersed boundary is widely adopted to study particle-laden flows. The effective hydrodynamic diameter of the particle is much or less overestimated by the original immersed boundary method, depending on the particle Reynolds number and grid resolution. In this paper, we propose an improved method to dynamically correct the effective hydrodynamic diameter by retracting inward the Lagrangian points to varying distances. The retraction distance is determined by querying a function fitted in this paper. The improved method is tested and validated by several cases, including falling of a spherical particle under gravity, the uniform flow past two stationary particles and the drafting-kissing-tumbling phenomenon of two settling particles. It turns out the improved method not only provides better results for the drag force but also predicts the flow field more accurately. What's more, due to the insensitivity to grid resolution, the improved method is suitable for simulating large-scale fluid–particle systems such as fluidized bed, which are computationally expensive.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): François P. Hamon, Martin Schreiber, Michael L. Minion〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉Efficient time integration schemes are necessary to capture the complex processes involved in atmospheric flows over long periods of time. In this work, we propose a high-order, implicit–explicit numerical scheme that combines Multi-Level Spectral Deferred Corrections (MLSDC) and the Spherical Harmonics (SH) transform to solve the wave-propagation problems arising from the shallow-water equations on the rotating sphere.〈/p〉 〈p〉The iterative temporal integration is based on a sequence of corrections distributed on coupled space–time levels to perform a significant portion of the calculations on a coarse representation of the problem and hence to reduce the time-to-solution while preserving accuracy. In our scheme, referred to as MLSDC-SH, the spatial discretization plays a key role in the efficiency of MLSDC, since the SH basis allows for consistent transfer functions between space–time levels that preserve important physical properties of the solution.〈/p〉 〈p〉We study the performance of the MLSDC-SH scheme with shallow-water test cases commonly used in numerical atmospheric modeling. We use this suite of test cases, which gradually adds more complexity to the nonlinear system of governing partial differential equations, to perform a detailed analysis of the accuracy of MLSDC-SH upon refinement in time. We illustrate the stability properties of MLSDC-SH and show that the proposed scheme achieves up to eighth-order convergence in time. Finally, we study the conditions in which MLSDC-SH achieves its theoretical speedup, and we show that it can significantly reduce the computational cost compared to single-level Spectral Deferred Corrections (SDC).〈/p〉 〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Lincheng Xu, Fang-Bao Tian, John Young, Joseph C.S. Lai〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉A novel computational framework which combines the lattice Boltzmann method (LBM) and an improved immersed boundary method (IBM) based on a dynamic geometry-adaptive Cartesian grid system is introduced for the fluid–structure interaction (FSI) problems at moderate and high Reynolds numbers. In this framework, the fluid dynamics is obtained by solving the discrete lattice Boltzmann equation. The boundary conditions at the fluid–structure interfaces are handled by an improved IBM based on a feedback scheme, which drives the predicted flow velocity (calculated after the LBM stream process without the IBM body force) near the immersed boundaries to match the solid velocity. In the present IBM, the feedback coefficient is mathematically derived and explicitly approximated. The Lagrangian force density is divided into two parts: one is the traction caused by the predicted flow velocity, and the other is caused by the acceleration of the immersed boundary. Such treatment significantly enhances the numerical stability for modelling FSI problems involving small structure-to-fluid mass ratios. A novel dynamic geometry-adaptive refinement is applied to provide fine resolution around the immersed geometries and coarse resolution in the far field. The overlapping grids between two adjacent refinements consist of two layers. In order to enhance the numerical stability, two-layer “ghost nodes” are generated within the immersed body domain which is a non-fluid area. The movement of fluid–structure interfaces only causes adding or removing grids at the boundaries of refinements and consequently a high mesh-update efficiency is guaranteed. Finally, large eddy simulation models are incorporated into the framework to model turbulent flows at relatively high Reynolds numbers. Several validation cases, including an impulsively started flow over a vertical plate, flow over stationary and oscillating cylinders, flow over flapping foils, flexible filaments in a uniform flow, turbulent flow over a wavy boundary, flow over a stationary sphere and a dragonfly in hovering flight, are conducted to verify the accuracy and fidelity of the present solver over a range of Reynolds numbers.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Carlos Pérez-Arancibia, Luiz M. Faria, Catalin Turc〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present an effective harmonic density interpolation method for the numerical evaluation of singular and nearly singular Laplace boundary integral operators and layer potentials in two and three spatial dimensions. The method relies on the use of Green's third identity and local Taylor-like interpolations of density functions in terms of harmonic polynomials. The proposed technique effectively regularizes the singularities present in boundary integral operators and layer potentials, and recasts the latter in terms of integrands that are bounded or even more regular, depending on the order of the density interpolation. The resulting boundary integrals can then be easily, accurately, and inexpensively evaluated by means of standard quadrature rules. A variety of numerical examples demonstrate the effectiveness of the technique when used in conjunction with the classical trapezoidal rule (to integrate over smooth curves) in two-dimensions, and with a Chebyshev-type quadrature rule (to integrate over surfaces given as unions of non-overlapping quadrilateral patches) in three-dimensions.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Kyle Gerard Felker, James M. Stone〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present a fourth-order accurate finite volume method for the solution of ideal magnetohydrodynamics (MHD). The numerical method combines high-order quadrature rules in the solution of semi-discrete formulations of hyperbolic conservation laws with the upwind constrained transport (UCT) framework to ensure that the divergence-free constraint of the magnetic field is satisfied. A novel implementation of UCT that uses the piecewise parabolic method (PPM) for the reconstruction of magnetic fields at cell corners in 2D is introduced. The resulting scheme can be expressed as the extension of the second-order accurate constrained transport (CT) Godunov-type scheme that is currently used in the Athena astrophysics code. After validating the base algorithm on a series of hydrodynamics test problems, we present the results of multidimensional MHD test problems which demonstrate formal fourth-order convergence for smooth problems, robustness for discontinuous problems, and improved accuracy relative to the second-order scheme.〈/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. Stanier, L. Chacón, G. Chen〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The quasi-neutral hybrid model with kinetic ions and fluid electrons is a promising approach for bridging the inherent multi-scale nature of many problems in space and laboratory plasmas. Here, a novel, implicit, particle-in-cell based scheme for the hybrid model is derived for fully 3D electromagnetic problems with multiple ion species, which features global mass, momentum and energy conservation. The scheme includes sub-cycling and orbit-averaging for the ions, with cell-centered finite differences and implicit midpoint time advance. To reduce discrete particle noise, the scheme allows arbitrary-order shape functions for the particle-mesh interpolations and the application of conservative binomial smoothing. The algorithm is verified for a number of test problems to demonstrate the correctness of the implementation, the unique conservation properties, and the favorable stability properties of the new scheme. In particular, there is no indication of unstable growth of the finite-grid instability for a population of cold ions drifting through a uniform spatial mesh, in a set-up where several commonly used non-conservative schemes are highly unstable.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): R. Chakir, Y. Maday, P. Parnaudeau〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Computation Fluid Dynamics (CFD) simulation has become a routine design tool for i) predicting accurately the thermal performances of electronics set ups and devices such as cooling system and ii) optimizing configurations. Although CFD simulations using discretization methods such as finite volume or finite element can be performed at different scales, from component/board levels to larger system, these classical discretization techniques can prove to be too costly and time consuming, especially in the case of optimization purposes where similar systems, with different design parameters have to be solved sequentially. The design parameters can be of geometric nature or related to the boundary conditions. This motivates our interest on model reduction and particularly on reduced basis methods. As is well documented in the literature, the offline/online implementation of the standard RB method (a Galerkin approach within the reduced basis space) requires to modify the original CFD calculation code, which for a commercial one may be problematic even impossible. For this reason, we have proposed in a previous paper, with an application to a simple scalar convection diffusion problem, an alternative non-intrusive reduced basis approach (NIRB) based on a two-grid finite element discretization. Here also the process is two stages: 〈em〉offline〈/em〉, the construction of the reduced basis is performed on a fine mesh; 〈em〉online〈/em〉 a new configuration is simulated using a coarse mesh. While such a coarse solution, can be computed quickly enough to be used in a rapid decision process, it is generally not accurate enough for practical use. In order to retrieve accuracy, we first project every such coarse solution into the reduced space, and then further improve them via a rectification technique. The purpose of this paper is to generalize the approach to a CFD configuration.〈/p〉〈/div〉

Description: 〈p〉Publication date: 1 January 2019〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 376〈/p〉 〈p〉Author(s): Eric J. Ching, Yu Lv, Peter Gnoffo, Michael Barnhardt, Matthias Ihme〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉This study is concerned with predicting surface heat transfer in viscous hypersonic flows using high-order discontinuous Galerkin (DG) methods. Currently, finite-volume (FV) schemes are most commonly employed for computing flows in which surface heat transfer is a target quantity; however, these schemes suffer from large sensitivities to a variety of factors, such as the inviscid flux function and the computational mesh. High-order DG methods offer advantages that can mitigate these sensitivities. As such, a simple and robust shock capturing method is developed for DG schemes. The method combines intraelement variations for shock detection with smooth artificial viscosity (AV) for shock stabilization. A parametric study is performed to evaluate the effects of AV on the solution. The shock capturing method is employed to accurately compute double Mach reflection and viscous hypersonic flows over a circular half-cylinder and a double cone, the latter of which involves a complex flow topology with multiple shock interactions and flow separation. Results show this methodology to be significantly less sensitive than FV schemes to mesh topology and inviscid flux function. Furthermore, quantitative comparisons with state-of-the-art FV calculations from an error vs. cost perspective are provided.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Gurpreet Singh, Gergina Pencheva, Mary F. Wheeler〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We present an approximate Jacobian approach for solving nonlinear, multiphase flow and transport problems in porous media. A backward Euler time discretization scheme is used prior to spatial discretization with a lowest order mixed finite element method (MFEM). This results in a fully implicit nonlinear algebraic system of equations. Conventionally, an exact Jacobian construction is employed during the Newton linearization to obtain a linear system of equations after spatial and temporal discretization. This fully coupled, monolithic linear system, usually in pressure and saturation (or concentration) unknowns, requires specialized preconditioners such as constrained pressure residual (CPR) or two stage preconditioner. These preconditioners operate on the linear system to decouple pressure and saturation (or concentration) degrees of freedom (DOF) in order to use existing linear solvers for positive definite (PD) matrices such as GMRES or AMG, to name a few. In this work, we present an alternative to two-stage preconditioning (or CPR) for solving the aforementioned monolithic system after Newton linearization. This approach relies upon a decoupling approximation for the pressure-saturation (or concentration) block sub-matrices, during Newton linearization, to obtain block diagonal sub-matrices. The resulting linear system is easily reduced, trivially eliminating these diagonal sub-matrices, to obtain a system in pressure DOF circumventing the need for specialized preconditioners. Further, the linear system has lesser DOF owing to the elimination of saturation (or concentration) unknowns. This nonlinear solver is demonstrated to be as accurate as the exact Jacobian approach, measured in terms of convergence of nonlinear residual to a desired tolerance for both methods. Our numerical results indicate a consistent computational speedup by a factor of approximately 1.32 to 4.0 for the two-phase flow model formulation under consideration. This is related to the DOF of the linear systems for the approximate and exact Jacobian approaches. For multicomponent flow and transport, this speedup is expected to be directly proportional to the number of concentration degrees of freedom. A number of field scale numerical simulations are also presented to demonstrate the efficacy of this approach for realistic problems.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Luis L. Bonilla, Ana Carpio, Manuel Carretero, Gema Duro, Mihaela Negreanu, Filippo Terragni〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We study a robust finite difference scheme for integrodifferential kinetic systems of Fokker–Planck type modeling tumor driven blood vessel growth. The scheme is of order one and enjoys positivity features. We analyze stability and convergence properties, and show that soliton-like asymptotic solutions are correctly captured. We also find good agreement with the solution of the original stochastic model from which the deterministic kinetic equations are derived working with ensemble averages. A numerical study clarifies the influence of velocity cut-offs on the solutions for exponentially decaying data.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Rakesh Kumar, Praveen Chandrashekar〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In the present work, we propose two new variants of fifth order finite difference WENO schemes of adaptive order. We compare our proposed schemes with other variants of WENO schemes with special emphasize on WENO-AO(5,3) scheme (Balsara et al., 2016) [3]. The first algorithm (WENO-AON(5,3)), involves the construction of a new simple smoothness indicator which reduces the computational cost of WENO-AO(5,3) scheme. Numerical experiments show that accuracy of WENO-AON(5,3) scheme is comparable to that of WENO-AO(5,3) scheme and resolution of solutions involving shock or other discontinuities is comparable to that of WENO-AO(5,3) scheme. The second algorithm denoted as WENO-AO(5,4,3), involves the inclusion of an extra cubic polynomial reconstruction in the base WENO-AO(5,3) scheme, which leads to a more accurate scheme. Extensive numerical experiments in 1D and 2D are performed, which shows that WENO-AO(5,4,3) scheme has better resolution near shocks or discontinuities among the considered WENO schemes with negligible increase in computational cost.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Van-Dang Nguyen, Johan Jansson, Johan Hoffman, Jing-Rebecca Li〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉 〈p〉The Bloch–Torrey equation describes the evolution of the spin (usually water proton) magnetization under the influence of applied magnetic field gradients and is commonly used in numerical simulations for diffusion MRI and NMR. Microscopic heterogeneity inside the imaging voxel is modeled by interfaces inside the simulation domain, where a discontinuity in the magnetization across the interfaces is produced via a permeability coefficient on the interfaces. To avoid having to simulate on a computational domain that is the size of an entire imaging voxel, which is often much larger than the scale of the microscopic heterogeneity as well as the mean spin diffusion displacement, smaller representative volumes of the imaging medium can be used as the simulation domain. In this case, the exterior boundaries of a representative volume either must be far away from the initial positions of the spins or suitable boundary conditions must be found to allow the movement of spins across these exterior boundaries.〈/p〉 〈p〉Many approaches have been taken to solve the Bloch–Torrey equation but an efficient high performance computing framework is still missing. In this paper, we present formulations of the interface as well as the exterior boundary conditions that are computationally efficient and suitable for arbitrary order finite elements and parallelization. In particular, the formulations are based on the partition of unity concept which allows for a discontinuous solution across interfaces conforming with the mesh with weak enforcement of real (in the case of interior interfaces) and artificial (in the case of exterior boundaries) permeability conditions as well as an operator splitting for the exterior boundary conditions. The method is straightforward to implement and it is available in FEniCS for moderate-scale simulations and in FEniCS-HPC for large-scale simulations. The order of accuracy of the resulting method is validated in numerical tests and a good scalability is shown for the parallel implementation. We show that the simulated dMRI signals offer good approximations to reference signals in cases where the latter are available and we performed simulations for a realistic model of a neuron to show that the method can be used for complex geometries.〈/p〉 〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Yongle Du〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉Co-existence of the physical and numerical boundary conditions makes implicit boundary treatment a particularly difficult problem in modern CFD simulations. Previous studies adopted space–time extrapolation or specially designed partial differential equations on the boundaries that are different from those of interior points. They are often formulated in terms of primitive variables, and are very challenging for complicated boundary types to be converted to, and implemented in, the conservative variables that are preferred in numerical simulations. More importantly, different boundary equations or different extrapolation techniques may compromise the stability, accuracy, or convergence rate of the A-stable schemes that are developed for interior points. A new methodology for implicit boundary treatment is proposed in this study. By introducing a simple correction matrix 〈em〉T〈/em〉, a set of generalized equations 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mo〉∂〈/mo〉〈mi〉Q〈/mi〉〈mo stretchy="false"〉/〈/mo〉〈mo〉∂〈/mo〉〈mi〉t〈/mi〉〈mo〉=〈/mo〉〈mo stretchy="false"〉(〈/mo〉〈mi〉I〈/mi〉〈mo〉+〈/mo〉〈mi〉T〈/mi〉〈mo stretchy="false"〉)〈/mo〉〈mi〉R〈/mi〉〈/math〉 are developed in terms of conservative variables. It is applicable for both the interior domain (〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mi〉T〈/mi〉〈mo〉=〈/mo〉〈mn〉0〈/mn〉〈/math〉) and the boundaries (〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si3.gif" overflow="scroll"〉〈mi〉T〈/mi〉〈mo〉≠〈/mo〉〈mn〉0〈/mn〉〈/math〉). It is in a partial differential equation form, satisfies the boundary conditions accurately, independent of the time and spatial discretizations. Any one-sided schemes can be used on the boundaries but still maintain the upwind property. Implicit solution techniques are made significantly easy to implement using, for example, the data-parallel lower–upper relation method and the Newton method (combined with the GMRES method for subsidiary iterations). Numerical experiments show that the proposed methodology produces stable simulations for very large CFL numbers and preserve the imposed boundary values accurately.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): John E. Ortiz G., Axelle Pillain, Lyes Rahmouni, Francesco P. Andriulli〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The symmetric formulation of the electroencephalography (EEG) forward problem is a well-known and widespread equation thanks to the high level of accuracy that it delivers. However, this equation is first kind in nature and gives rise to ill-conditioned problems when the discretization density or the brain conductivity contrast increases, resulting in numerical instabilities and increasingly slow solutions. This work addresses and solves this problem by proposing a new regularized symmetric formulation. The new scheme is obtained by leveraging on Calderon identities which allow to introduce a dual symmetric equation that, combined with the standard one, results in a second kind operator which is both stable and well-conditioned under all the above mentioned conditions. The new formulation presented here can be easily integrated into existing EEG imaging packages since it can be obtained with the same computational technology required by the standard symmetric formulation. The performance of the new scheme is substantiated by both theoretical developments and numerical results which corroborate the theory and show the practical impact of the new technique.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): Arnout M.P. Boelens, Daniele Venturi, Daniel M. Tartakovsky〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉High-dimensional partial-differential equations (PDEs) arise in a number of fields of science and engineering, where they are used to describe the evolution of joint probability functions. Their examples include the Boltzmann and Fokker–Planck equations. We develop new parallel algorithms to solve such high-dimensional PDEs. The algorithms are based on canonical and hierarchical numerical tensor methods combined with alternating least squares and hierarchical singular value decomposition. Both implicit and explicit integration schemes are presented and discussed. We demonstrate the accuracy and efficiency of the proposed new algorithms in computing the numerical solution to both an advection equation in six variables plus time and a linearized version of the Boltzmann equation.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): J. Shipton, T.H. Gibson, C.J. Cotter〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉We describe a compatible finite element discretisation for the shallow water equations on the rotating sphere, concentrating on integrating consistent upwind stabilisation into the framework. Although the prognostic variables are velocity and layer depth, the discretisation has a diagnostic potential vorticity that satisfies a stable upwinded advection equation through a Taylor–Galerkin scheme; this provides a mechanism for dissipating enstrophy at the gridscale whilst retaining optimal order consistency. We also use upwind discontinuous Galerkin schemes for the transport of layer depth. These transport schemes are incorporated into a semi-implicit formulation that is facilitated by a hybridisation method for solving the resulting mixed Helmholtz equation. We demonstrate that our discretisation achieves the expected second order convergence and provide results from some standard rotating sphere test problems.〈/p〉〈/div〉

Description: 〈p〉Publication date: 15 December 2018〈/p〉 〈p〉〈b〉Source:〈/b〉 Journal of Computational Physics, Volume 375〈/p〉 〈p〉Author(s): T.A. Biala, A.Q.M. Khaliq〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉In this paper, we develop a time stepping scheme for solving nonlinear time–space fractional partial differential equations (PDEs). In space, we use the matrix transfer technique to discretize the PDEs and obtain a system of nonlinear time-fractional differential equations. The developed scheme is similar to the Crank–Nicholson scheme for integer order PDEs and are shown to be of order 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mn〉1〈/mn〉〈mo〉+〈/mo〉〈mi〉α〈/mi〉〈/math〉 in time where 〈em〉α〈/em〉 is the order of the time derivative described in the Caputo sense. The solution of the PDE at any point 〈em〉i〈/em〉 in the 〈em〉t〈/em〉-stencil depends not only on the solution at the point 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mi〉i〈/mi〉〈mo〉−〈/mo〉〈mn〉1〈/mn〉〈/math〉 but on all previous solutions (memory effect). Thus, the implementation of schemes for such fractional PDEs, for long time interval, can be time consuming. This is basically due to the computation and re-computation of the history term at each time step. We lessen this computational time by implementing three parallel versions of the algorithm. The shared memory systems (OpenMP) and the distributed memory systems (MPI) are used for implementing the parallel algorithms. A third parallel version uses both the shared and distributed memory systems (Hybrid version). The advantages of the parallel algorithms over the sequential algorithm are discussed. The merits and demerits of each parallel versions of the algorithm over the others are examined. Numerical simulations are performed to support our theoretical observations.〈/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. Esmaily, J.A.K. Horwitz〈/p〉 〈h5〉Abstract〈/h5〉 〈div〉〈p〉The accuracy of Lagrangian point-particle models for simulation of particle-laden flows may degrade when the particle and fluid momentum equations are two-way coupled. The exchange of force between the fluid and particle changes the velocity of the fluid at the location of the particle, thereby modifying the slip velocity and producing an erroneous prediction of coupling forces between fluid and particle. In this article, we propose a correction scheme to reduce this error and predict the undisturbed fluid velocity accurately. Conceptually, in this method, the computation cell is treated as a solid object immersed in the fluid that is subjected to the two-way coupling force and dragged at a velocity that is identical to the disturbance created by the particle. The proposed scheme is generic as it can be applied to unstructured grids with arbitrary geometry, particles that have different size and density, and arbitrary interpolation scheme. In its crudest form for isotropic grids, the present correction scheme reduces to dividing the Stokes drag by 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si1.gif" overflow="scroll"〉〈mn〉1〈/mn〉〈mo〉−〈/mo〉〈mn〉0.75〈/mn〉〈mi mathvariant="normal"〉Λ〈/mi〉〈/math〉 for 〈math xmlns:mml="http://www.w3.org/1998/Math/MathML" altimg="si2.gif" overflow="scroll"〉〈mi mathvariant="normal"〉Λ〈/mi〉〈mo〉≤〈/mo〉〈mn〉1〈/mn〉〈/math〉, where Λ is the ratio of the particle diameter to the grid size. The accuracy of the proposed scheme is evaluated by comparing the computed settling velocity of an individual and pair of particles under gravity on anisotropic rectilinear grids against analytical solutions. This comparison shows up to two orders of magnitude reduction in error in cases where the particle is up to 5 times larger than the grid that may have an aspect ratio of over 10. Furthermore, a comparison against a particle-resolved simulation of decaying isotropic turbulence demonstrates the excellent accuracy of the proposed scheme.〈/p〉〈/div〉