Grid Refinement Research Articles

A 2D-SWEs framework for efficient catchment-scale simulations: Hydrodynamic scaling properties of river networks and implications for non-uniform grids generation

The application of two-dimensional shallow-water equations models (2D-SWEs) for the description of hydrodynamic-based surface runoff computations is becoming a reference approach in rainfall-runoff simulations at the catchment scale. Due to their ability in generation of flow patterns throughout the basin, they can be used not only as an advanced method for flood mapping studies and hazard assessment but also as an innovative tool for the analysis of river drainage networks, opening new perspectives for several environmental processes. In particular, in this work we put the river networks in a 2D-SWEs framework, meaning that the traditional tree-like fluvial structure, represented by a skeleton composed of a set of lines, is replaced by a collection of points discretizing the 2-D geometry of the river structure itself, for which the values of the hydrodynamic values are provided by the numerical simulations. This approach is used here to derive a new scaling property that relates the specific discharge threshold, used to identify the river network cells, to the total areas of the network cells themselves. The hydrodynamic and geomorphological interpretation of this power law function and the influence of grid resolution, on some relevant parameters of this curve, have inspired the development of a heuristic procedure for non-uniform grid generation, able to detect the most hydrodynamically active areas of the basins for which the grid refinement process makes sense. Moreover, information related to how much grid refinement is needed is provided as well. The performances of this procedure are very promising in terms of accuracy of simulated discharges, hydrodynamic behaviour of the river network and flooded areas, reducing significantly the computational times in respect to the use of fine uniform grids.

Robust DOA estimation for burst impulsive noise

This paper describes a robust sparse Bayesian learning approach for DOA estimation in the presence of impulsive noise. In practical situations, impulses in the noise often appear in bursts, resulting in structured rather than independent sparsity, which can be exploited to enhance DOA estimation performance. However, there is no method available for capturing the realistic burst structure of impulsive noise in the literature. In this paper, we model the structured impulsive noise with a novel two-dimensional burst prior and develop a robust variational Bayesian inference (VBI) framework for DOA estimation under burst impulsive noise, and an improved grid refinement method is further proposed to adjust grid points of the dictionary matrix, which can combat the off-grid gap more efficiently. The superior impulse-resistant DOA estimation performance and computational efficiency of our proposed method are verified by simulation results.

Plasma edge simulations including realistic wall geometry with SOLPS-ITER

Abstract In plasma edge simulations using the SOLPS-ITER code, the simulated Scrape-Off Layer plasma domain has historically been restricted to magnetic flux surfaces contacting divertor targets at both ends. We present here a newly developed numerical solver for the B2.5 plasma solver in SOLPS-ITER, allowing the numerical grid to be extended to the true vessel boundaries. The new, unstructured Finite Volume scheme can deal with arbitrary grids and magnetic topologies in the 2D poloidal plane. It includes a correct numerical treatment of possibly misaligned faces and cells w.r.t. the magnetic field to cope with, for example, strong divertor target shaping. The solver combines the benefits of an accurate numerical separation of fast parallel and slow radial transport, with a realistic description of the wall geometry, and the possibility of local grid refinement to capture sharp features in the Scrape-Off Layer flows. Generalized sheath boundary conditions are presented that can be imposed at all vessel boundaries, removing an important modeling uncertainty related to the specification of ad hoc decay length boundary conditions at the outer flux surfaces. The resulting model is applied to an AUG single-null case, a standard benchmark case for SOLPS-ITER. We analyze in particular the impact of the extended plasma model on upstream and divertor plasma conditions, and the improved predictions of heat and particle loads to the main chamber wall. The extended solver also allows for a much improved qualitative agreement between fluid and kinetic neutral simulations, because the fluid neutral solution, which is obtained on the plasma grid, now also extends to the true main chamber and divertor vessel boundaries.

Kinetic Functions for Nonclassical Shocks, Entropy Stability, and Discrete Summation by Parts

We study nonlinear hyperbolic conservation laws with non-convex flux in one space dimension and, for a broad class of numerical methods based on summation by parts operators, we compute numerically the kinetic functions associated with each scheme. As established by LeFloch and collaborators, kinetic functions (for continuous or discrete models) uniquely characterize the macro-scale dynamics of small-scale dependent, undercompressive, nonclassical shock waves. We show here that various entropy-dissipative numerical schemes can yield nonclassical solutions containing classical shocks, including Fourier methods with (super-) spectral viscosity, finite difference schemes with artificial dissipation, discontinuous Galerkin schemes with or without modal filtering, and TeCNO schemes. We demonstrate numerically that entropy stability does not imply uniqueness of the limiting numerical solutions for scalar conservation laws in one space dimension, and we compute the associated kinetic functions in order to distinguish between these schemes. In addition, we design entropy-dissipative schemes for the Keyfitz-Kranzer system whose solutions are measures with delta shocks. This system illustrates the fact that entropy stability does not imply boundedness under grid refinement.

An analytical tool to forecast horizontal well injectivity in viscous oil polymer floods

Injection rate is a key factor for successful polymer floods of viscous oil reservoirs. Injectivity during polymer flooding is lower compared to waterflooding due to the viscosity of polymer solution. Accurate estimation of injection rate/pressure build-up due to polymer injection is essential to evaluation and design of polymer flood projects. Typically, a short-term polymer injectivity test is performed on a single well, but it is hard to extrapolate future behavior of polymer injectivity because the tests are often too short. Reservoir simulation can be used to forecast injectivity, but the grid size and radial flow often affect the calculated injectivity. Analytical methods provide simple and quick ways to estimate injectivity using realistic reservoir conditions. The goal of this work is to develop an analytical tool and demonstrate its applicability by comparison with field data. An analytical solution is developed for the injectivity of horizontal wells in polymer floods in viscous oil reservoirs. The analytical solution is validated by comparing it to a numerical field-scale box model simulation. The pressure profiles are calculated using spreadsheets and compared against the data from several fields with increasing degrees of complexity. In each case, the analytical-spreadsheet method provides close matches of injection pressure data. This analytical tool can be used for preliminary design of polymer floods before getting into more complex and time-consuming reservoir simulations. • Developed an analytical solution for injectivity in viscous oil polymer floods. • Compared the analytical solution to numerical simulation with grid refinement. • Compared analytical solution to field data from 5 wells. • Novel tool to design and optimize polymer floods in viscous oil reservoirs.

Performance Measurements Within Asynchronous Task-Based Runtime Systems: A Double White Dwarf Merger as an Application

Analyzing performance within asynchronous many-task-based runtime systems is challenging because millions of tasks are launched concurrently. Especially for long-term runs the amount of data collected becomes overwhelming. We study HPX and its performance-counter framework and APEX to collect performance data and energy consumption. We added HPX application-specific performance counters to the Octo-Tiger full 3D AMR astrophysics application. This enables the combined visualization of physical and performance data to highlight bottlenecks with respect to different solvers. We examine the overhead introduced by these measurements, which is around 1%, with respect to the overall application runtime. We perform a convergence study for four different levels of refinement and analyze the application's performance with respect to adaptive grid refinement. The measurements' overheads are small, enabling the combined use of performance data and physical properties with the goal of improving the code's performance. All of these measurements were obtained on NERSC's Cori, Louisiana Optical Network Infrastructure's QueenBee2, and Indiana University's Big Red 3.

Implementation of a hybrid Lagrangian filtered density function–large eddy simulation methodology in a dynamic adaptive mesh refinement environment

Multispecies mixing processes play an important role in many engineering, biological, and environmental applications. Since simulating mixing flows can be useful to understand its physics and to study industrial issues, this work aims to develop the basis of a methodology able to simulate the physics of multiple-species mixing flows, using a hybrid large eddy simulation/Lagrangian filtered density function (FDF) method on an adaptive, block-structured mesh. A computational model of notional particles transport on a distributed processing environment is built using a parallel Lagrangian map. This map connects the Lagrangian information with the Eulerian framework of the in-house code MFSim, in which transport equations are solved. The Lagrangian composition FDF method, through the Monte Carlo technique, performs algebraic calculations over an ensemble of notional particles and provides composition fields statistically equivalent to those obtained by finite volume numerical solution of partial differential equations. Finally, to maintain high accuracy in the system of stochastic differential equations solver when an adaptive mesh refinement environment is used, a methodology for ensuring mass conservation is developed to preserve at least the statistical moments up to order two, even in the case of annihilation or cloning of a large number of notional particles in one time step, ensuring the applicability of Lagrangian FDF methods in dynamically adaptive grid refinement.

Transition region adaptive conduction (TRAC) in multidimensional magnetohydrodynamic simulations

Context. In solar physics, a severe numerical challenge for modern simulations is properly representing a transition region between the million-degree hot corona and a much cooler plasma of about 10 000 K (e.g., the upper chromosphere or a prominence). In previous 1D hydrodynamic simulations, the transition region adaptive conduction (TRAC) method has been proven to capture aspects better that are related to mass evaporation and energy exchange. Aims. We aim to extend this method to fully multidimensional magnetohydrodynamic (MHD) settings, as required for any realistic application in the solar atmosphere. Because modern MHD simulation tools efficiently exploit parallel supercomputers and can handle automated grid refinement, we design strategies for any-dimensional block grid-adaptive MHD simulations. Methods. We propose two different strategies and demonstrate their working with our open-source MPI-AMRVAC code. We benchmark both strategies on 2D prominence formation based on the evaporation–condensation scenario, where chromospheric plasma is evaporated through the transition region and then is collected and ultimately condenses in the corona. Results. A field-line-based TRACL method and a block-based TRACB method are introduced and compared in block grid-adaptive 2D MHD simulations. Both methods yield similar results and are shown to satisfactorily correct the underestimated chromospheric evaporation, which comes from a poor spatial resolution in the transition region. Conclusions. Because fully resolving the transition region in multidimensional MHD settings is virtually impossible, TRACB or TRACL methods will be needed in any 2D or 3D simulations involving transition region physics.

USM3D Analyses of Transonic Boundary Layer Thickener Configurations Applied to Fuselage

A wind tunnel test was performed in the 0.3-Meter Cryogenic Transonic Tunnel at the NASA Langley Research Center to investigate the impact of various boundary layer thickeners on the height and shape of the boundary layer on the aft portion of a half-fuselage model. The present paper discusses USM3D analyses that were performed in support of the wind tunnel test. The results consisted of comparisons between USM3D and experimental data for the baseline configuration and two boundary layer thickening configurations. The comparisons considered surface pressure distributions, boundary layer profiles, and calculated boundary layer thicknesses. Additionally, a grid refinement study was performed for each of the three configurations, which illustrated general improvement in the comparisons with increasing grid refinement. Overall, the comparisons show favorable agreement between the USM3D predictions and the experimental data.

Performance Analysis of FEM Solvers on Practical Electromagnetic Problems

The paper presents a comparative analysis of different commercial and academic software. The comparison aims to examine how the integrated adaptive grid refinement methodologies can deal with challenging, electromagnetic-field related problems. For this comparison, two bench-mark problems were examined in the paper. The first example is a solution of an L-shape domain like test problem, which has a singularity at a certain point in the geometry. The second problem is an induction heated aluminum rod, which accurate solution needs to solve non-linear, coupled physical fields. The accurate solution of this problem requires applying adaptive mesh generation strategies or applying a very fine mesh in the electromagnetic domain, which can significantly increase the computational complexity. The results show that the fully-hp adaptive meshing strategies, which are integrated into Agros Suite, can significantly reduce the task's computational complexity compared to the automatic h-adaptivity, which is part of the examined, popular commercial solvers.

Computational methods for pore-scale simulation of rarefied gas flow

Direct simulation at the pore-scale is crucial to unravel rarefaction effect on gas transport in tight porous media. To satisfy the dual demands on modeling accuracy and computation effort, an appropriate method must be chosen. This work, therefore, evaluates four numerical methods for pore-scale rarefied gas flows in a two-dimensional (2D) model porous media over a wide range of rarefaction. These methods are the incompressible Navier-Stokes equation with the first-order velocity-slip boundary condition, and three gas-kinetic solvers i.e. a finite-difference (FD) iterative solver for the linearized Bhatnagar, Gross and Crook (BGK) model kinetic equation, a finite-volume (FV) solver for the non-linear Shakhov model, and an open-source direct simulation Monte Carlos (DSMC) solver. The benchmark cases cover the characteristic Knudsen number ranging from 0.0231 to 4.62 while the characteristic Reynolds numbers are kept to be less than 1.0. All the solvers are developed in OpenFOAM, except the FD solver, allowing us to investigate the effect of local grid refinements using the automatic Cartesian grid generator in OpenFOAM. The flow fields and the apparent permeabilities predicted by all the solvers have been compared in detail. Besides, the computational time of these solvers is measured and analyzed to demonstrate the relative cost of the three kinetic solvers. It is found that the FD solver is the most efficient one and gives accurate results over the whole range of Knudsen number. Finally, this study also evaluates the feasibility of the recently developed discrete unified gas-kinetic scheme (DUGKS), the algorithm in the aforementioned non-linear Shakhov model equation solver, for pore-scale rarefied gas flow. It is found the results predicted by this solver agree well with the other two kinetic solvers.

Advance in compact structure‐preserving manner to the Rosenau–Kawahara model of shallow‐water wave

In this paper, we pursue further analysis of the performance of a compact structure‐preserving finite difference scheme. The convergence, stability, and accuracy of the approximate solution with respect to grid refinement are discussed. The compact difference approach to precisely preserving invariants on any time‐space regions gives a three‐level linear‐implicit scheme with the spatial accuracy, found to be fourth order on a uniform grid. The method is verified by comparison with a solution of the Rosenau–Kawahara equation just obtained with second‐order finite difference schemes recently. Also, the efficiency of the present algorithm is confirmed by simulations of the problem at a long time. Details of CPU time are examined in order to assess the usefulness of the compact scheme for determining an approximate solution.

Taitoko, an advanced code for tsunami propagation, developed at the French Tsunami Warning Centers

This paper presents the new version of the operational code used at the Tsunami Warning Centers in French Polynesia and in mainland France. This nested multigrid code solves either the nonlinear shallow water or the weakly nonlinear Boussinesq equations in a spherical or Cartesian coordinates system. The objective is to simulate the propagation of tsunamis both in the context of tsunami early warning and in the context of research. The originality of this code is to simulate the tsunami propagation through a series of nested grids in which the numerical scheme is specified by the user. Three standard finite-difference methods and two finite volume schemes are available to take into account the dispersive nature of waves or to capture flow discontinuities. The performance of the five numerical schemes has been assessed through idealized 1D and 2D test cases. The convergence rate with grid refinement is evaluated for three of these tests. Diffusion and dispersion of the five schemes are evaluated through the 1D far-field propagation in a channel of constant depth. Finally, two-dimensional simulations of the 2011 Japanese tsunami are performed considering a homogeneous tectonic source. Coastal inundation is simulated in the bay of Kamaishi where run-up heights reached 10 m. The ability of the five schemes to propagate trans-Pacific tsunamis is evaluated through comparisons of time series between simulations and observations at DART buoys. In the context of tsunami early warning, numerical results show that standard finite-difference methods reproduce most of the observed time gauges and that numerical errors associated to these schemes are acceptable.

Evaluation of Length Scales and Meshing Requirements for Resolving Two-Phase Flow Regime Transitions Using the Level Set Method

Abstract New criteria for fully resolving two-phase flow regime transitions using direct numerical simulation with the level set method for interface capturing are proposed. A series of flows chosen to capture small scale interface phenomena are simulated at different grid refinements. These cases include droplet deformation and breakup in a simple shear field, the thin film around a Taylor bubble, and the rise of a bubble toward a free surface. These cases cover the major small scale phenomena observed in two-phase flows: internal recirculation, interface curvature, interface snapping, flow of liquid in thin films, and drainage/snapping of thin films. The results from these simulations and their associated grid studies were used to develop new meshing requirements for simulation of two-phase flow using interface capturing methods, in particular the level set method. When applicable, the code used in this work, PHASTA, was compared to experiments in order to contribute to the ongoing validation process of the code. Results show that when the solver meets these criteria, with the exception of resolving the nanometer scale liquid film between coalescing bubbles, the code is capable of accurately simulating interface topology changes.

Improved Calculation of Wellblock Pressures for Numerical Simulation of Non-Newtonian Polymer Injection

Summary The viability of any enhanced-oil-recovery project depends on the ability to inject the displacing fluid at an economic rate. This is typically evaluated using finite-volume numerical simulation. These simulators calculate injectivity using the Peaceman method (Peaceman 1978), which assumes that flow is Newtonian. Most polymer solutions exhibit some degree of non-Newtonian behavior resulting in a changing polymer viscosity with distance from the injection well. For shear-thinning polymer solutions, conventional simulations can overpredict injection-well bottomhole pressure (BHP) by several hundred psi, unless a computationally costly local grid refinement is used in the near-wellbore region. We show theoretically and numerically that the Peaceman pressure-equivalent radius, based on Darcy flow, is not correct when fluids are shear thinning, and derive an analytical expression for calculating the correct radius. The expression does not depend on any particular functional relationship between polymer-solution viscosity and velocity. We test it using the relationship described by the Meter equation (Meter and Bird 1964) and the Cannella et al. (1988) correlation. Numerical tests indicate that the solution provides a significant improvement in the accuracy of BHP calculations for conventional numerical simulation, reducing or removing the need for expensive local grid refinement around the well when simulating the injection of fluids with shear-thinning non-Newtonian rheology.

High performance Wannier interpolation of Berry curvature and related quantities with WannierBerri code

Wannier interpolation is a powerful tool for performing Brillouin zone integrals over dense grids of k points, which are essential to evaluate such quantities as the intrinsic anomalous Hall conductivity or Boltzmann transport coefficients. However, more complex physical problems and materials create harder numerical challenges, and computations with the existing codes become very expensive, which often prevents reaching the desired accuracy. In this article, I present a series of methods that boost the speed of Wannier interpolation by several orders of magnitude. They include a combination of fast and slow Fourier transforms, explicit use of symmetries, and recursive adaptive grid refinement among others. The proposed methodology has been implemented in the python code WannierBerri, which also aims to serve as a convenient platform for the future development of interpolation schemes for other phenomena.

Unravelling cosmic velocity flows: a Helmholtz–Hodge decomposition algorithm for cosmological simulations

In the context of intra-cluster medium turbulence, it is essential to be able to split the turbulent velocity field in a compressive and a solenoidal component. We describe and implement a new method for this aim, i.e., performing a Helmholtz–Hodge decomposition, in multi-grid, multi-resolution descriptions, focusing on (but not being restricted to) the outputs of AMR cosmological simulations. The method is based on solving elliptic equations for a scalar and a vector potential, from which the compressive and the solenoidal velocity fields, respectively, are derived through differentiation. These equations are addressed using a combination of Fourier (for the base grid) and iterative (for the refinement grids) methods. We present several idealised tests for our implementation, reporting typical median errors in the order of 1‰-1%, and with 95-percentile errors below a few percents. Additionally, we also apply the code to the outcomes of a cosmological simulation, achieving similar accuracy at all resolutions, even in the case of highly non-linear velocity fields. We finally take a closer look to the decomposition of the velocity field around a massive galaxy cluster.

Dynamic effects of fault modeling on stair-step and corner-point grids

Fault modeling has become an integral element of reservoir simulation for structurally complex reservoirs. Modeling of faults in general has major implications for the simulation grid. In turn, the grid quality control is very important in order to attain accurate simulation results. We investigate the dynamic effects of using stair-step grid (SSG) and corner-point grid (CPG) approaches for fault modeling from the perspective of dynamic reservoir performance forecasting. We have performed a number of grid convergence and grid-type sensitivity studies for a variety of simple, yet intuitive faulted flow simulation problems with gradually increasing complexity. We have also explored the added value of the multipoint flux approximation (MPFA) method over the conventional two-point flux approximation (TPFA) to increase the accuracy of reservoir simulation results obtained on CPGs. Effects of fault seal modeling on grid-resolution convergence and grid-type sensitivity have also been briefly examined. For simple geometries, both SSG and CPG can be used for fault modeling with similar accuracy in conjunction with the pillar-grid approach. This is evidenced by the fact that simulation results from SSG and CPG converge to identical solutions. SSG and CPG yield different results for more complex geometries. Simulation results approach to a converged solution for relatively fine SSGs. However, a SSG only provides an approximation to the fault geometry and reservoir volumes when the grid is coarse. On the other hand, non-orthogonality errors are increasingly evident in relatively more complex faulted models on CPGs and such errors cannot be addressed by grid refinement. It has been observed that MPFA partially addresses the discretization errors on non-orthogonal grids but only from the total flux accuracy perspective. However, transport related errors are still evident. Grid convergence behaviors and grid effects are quite similar with or without fault seal modeling (i.e., dedicated fault-zone modeling by use of scaled-up seal factors) for simple geometries. However, in more complex test cases, we have observed that it is more difficult to achieve converged results in conjunction with fault seal modeling due to increased heterogeneity of the underlying problem.

Numerical simulation of boiling on unstructured grids

A numerical method for the simulation of multiphase flows with phase change on unstructured grids is presented. Based on the work of Tanguy et al. [31] using a Level Set/Ghost Fluid Method coupling for two-dimensional axisymmetric cartesian grids, we extend the method to two- and three-dimensional unstructured grids with the aim of taking a step towards realistic boiling simulations in industrial context. The mass transfer rate at the interface accounting for phase change is computed as a function of the liquid and vapor heat fluxes at the interface and of the latent heat of the fluid by means of a new framework improving the accuracy of differential operators. The mass transfer rate is then used in the projection method to solve Navier-Stokes equations and in the advection equation of the Level Set function to account for the interface movement due to phase change. An immersed Dirichlet boundary condition is imposed at the interface to ensure that boiling always occurs at saturation temperature. We demonstrate the accuracy of our method first in the case of a static growing bubble with a fixed constant mass transfer rate, and second with a mass transfer rate computed from the thermal fluxes at the interface. In both cases, the bubble radius at final time converges with grid refinement towards the theoretical value.

New approximations for one-dimensional 3-point and two-dimensional 5-point compact integrated RBF stencils

This paper presents some new compact approximation stencils based on integrated radial basis functions (IRBFs) for numerically solving second-order elliptic differential problems on Cartesian grids. Higher-order IRBF schemes are employed to approximate the field/dependent variable. The IRBF approximations in each direction are based on 3 points and constructed independently, where derivatives of the second, third, fourth, fifth and sixth orders along the grid line are enforced at the two end-points. The imposed nodal derivative values are simply acquired through a Picard-type iteration scheme. The stencil is made up of 3 points and 5 points for 1D and 2D discretisations, respectively. Numerical results show that the proposed stencils yield a high rate of convergence with respect to grid refinement, e.g. up to the 13th order for 1D problems and to the 9th order for 2D problems.