A Scheme for Generating Unstructured Grids on Spheres with Application to Parallel Computation

We present a comprehensive study of velocity interpolation methods in polygons. These methods are often used as postprocessing procedures for numerical schemes that do not directly calculate the velocity field but only provide cell boundary flux conditions, such as the finite volume schemes. These methods extend the widely used velocity interpolation algorithms, such as the Pollock's algorithm, to more complex geometries such as perpendicular bisection (PEBI) grids, unstructured triangular grids and grids with local refinement. Once the velocity field is interpolated, streamline trajectories and time of flight along the streamlines can be calculated for reservoir simulation, model calibration and waterflood management, for instance. These velocity interpolation methods assume known lower‐order or higher‐order cell boundary fluxes, which satisfy global mass conservation and normal flux continuity. However, they differ in the interpolation of velocities within the interior of the cells. The interpolating velocity may be locally conservative or nonconservative, continuous or discontinuous, lower order or higher order. Results show that the interpolated velocity field has to be locally conservative in order to guarantee the correct volumetric transformation for the calculated streamlines and the time of flight. Velocity continuity is not as important as local conservation for the purpose of streamline applications. Compared to higher‐order interpolation for the streamline trajectories, lower‐order interpolation has the advantage of an analytic solution and an efficient implementation. Based on our analysis, we recommend a lower‐order locally conservative method for the most robust and numerically efficient calculation of streamline trajectories on unstructured grids.

Effective workflows for translating earth models into simulation models require grids that preserve geologic accuracy, offer flexible resolution control, integrate tightly with upscaling, and can be generated easily. Corner-point grids and pillar-based unstructured grids fail to satisfy these objectives; hence, a truly 3D unstructured approach is required. This paper describes unstructured cut-cell gridding tools that address these needs and improve the integration of our overall reservoir modeling workflows. The construction of simulation grids begins with the geologic model: a numerical representation of the reservoir structure, stratigraphy, and properties. Our gridding utilizes a geochronological (GeoChron) map from physical coordinates to an unfaulted and unfolded depositional coordinate system. The mapping is represented implicitly on a tetrahedral mesh that conforms to faults, and it facilitates accurate geostatistical modeling of static depositional properties. In the simplest use case, we create an explicit representation of the geologic model as an unstructured polyhedral grid. Away from faults and other discontinuities, the cells are hexahedral, highly orthogonal, and arranged in a structured manner. Geometric cutting operations create general polyhedra adjacent to faults and explicit contact polygons across faults. The conversion of implicit models to explicit grids is conceptually straightforward, but the implementation is nontrivial due to the limitations of finite precision arithmetic and the need to remove small cells formed in the cutting process. In practice, simulation grids are often constructed at coarser resolutions than earth models. Our implementation of local grid coarsening and refinement exploits the flexibility of unstructured grids to minimize upscaling errors and preserve critical geologic features. Because the simulation grid and the geologic model are constructed using the same mapping, fine cells can be nested exactly inside coarse cells. Therefore, flow-based upscaling can be applied efficiently without resampling onto temporary local grids. This paper describes algorithms and data structures for constructing, storing, and simulating cut-cell grids. Examples illustrate accurate modeling of normal faults, y-faults, overturned layers, and complex stratigraphy. Flow results, including a field sector model, show the suitability of cut-cell grids for simulation.

A computational scheme for computing the water surface wave propagation, based by the nonlinear shallow water equations is described. The Taylor-Galerkin method is applied for discretization on time. Mixed Modified Finite Volume Method on unstructured triangular grid is used for spatial approximation. The linear basis functions are used for approximation wave surface, quadratic basis functions for approximation velocity fields for elimination non-physical spatial oscillation.

In reservoir simulation the pressure is the solution of an elliptic equation. It follows from the maximum principle that this equation satisfies a monotonicity property. This property should be preserved when discretising the pressure equation. If this is not the case, the pressure solution may have false internal oscillations and extrema on no-flow boundaries. Thus, there is a need for sufficient conditions for the discretization methods to be monotone. These conditions will depend on the permeability tensor and the grid. Previously monotonicity for control volume methods has been studied on quadrilateral grids and on hexagonal grids. The former was done for general methods which reproduces linear potential fields, while the latter was done in a Control Volume Finite Element setting. These analyses have given sharp sufficient conditions for the discretisation methods to be monotone. However, for methods whose cell stencils include cells which do not have any edges common with the central cell in the discretisation scheme, no work has been done on unstructured grids. In this work, we study monotonicity on triangular grids for control volume methods which are exact for linear potential fields. We derive sufficient conditions for monotonicity of the MPFA-O and -L methods. The found monotonicity regions for the MPFA methods are also tested numerically. The tests are done both on uniform grids in homogeneous media, and on perturbed grids, which corresponds to heterogeneous media. The investigations are done for single phase flow only. However, the results are relevant for multiphase simulations. The results obtained in this work may be utilised in grid generation. In this way we can construct grids where the discretisation of the pressure equation is guaranteed to be monotone.

The Parallel Block Jacobi (PBJ) spatial domain decomposition is well suited for implementation on massively parallel computers to solve the neutron transport equation on unstructured grids due to the simple scheduling policy that arises from the PBJ’s iterative asynchronicity. The Parallel Block Jacobi-Integral Transport Matrix Method (PBJ-ITMM) is an iterative method that utilizes the PBJ decomposition and resolves local within-group scattering in a single iteration, but requires a matrix-vector iterative solution. This work details the development, implementation, and testing of the novel Green’s Function ITMM Construction (GFIC) algorithm. The GFIC constructs the matrices required for the PBJ-ITMM’s iterative solution on unstructured grids, utilizing the physical interpretation of these matrices as discretized response functions to create a local problem with a Green’s Function–like source. Conducting a set of mesh sweeps over all angles on this local problem yields the ITMM matrix elements. On unstructured grids, this approach utilizes the kernel calculation and fundamental solution algorithm present in an existing transport code, thus avoiding reimplementation of code functionality. Using the GFIC, the PBJ-ITMM is implemented in THOR, a tetrahedral mesh transport code, along with the Inexact Parallel Block Jacobi (IPBJ) method for performance comparison. This comparison involves strong and weak scaling studies of the Godiva and C5G7 benchmark problems using up to 32 768 processors. These studies establish that the PBJ-ITMM executes faster than the IPBJ when the number of cells per subdomain falls below a problem-dependent threshold, ~128 cells for Godiva, >256 cells for C5G7. The largest problem tested, comprising more than 6.8 billion unknowns, solves in <30 min with the IPBJ and <20 min with the PBJ-ITMM, using 32 768 processors. These results demonstrate the PBJ-ITMM as a viable approach for solving neutron transport problems on unstructured grids using massively parallel computers. Additionally, this study illustrates the range of number of cells per subdomain over which this method is favorable.

Numerical analysis of electrical defibrillation: The parallel approach

A new code and methodology are introduced for solving the general relativistic magnetohydrodynamic (GRMHD) equations in fixed background spacetimes using time-explicit, finite-volume discretization. The code has options for solving the GRMHD equations using traditional artificial-viscosity (AV) or non-oscillatory central difference (NOCD) methods, or a new extended AV (eAV) scheme using artificial-viscosity together with a dual energy-flux-conserving formulation. The dual energy approach allows for accurate modeling of highly relativistic flows at boost factors well beyond what has been achieved to date by standard artificial viscosity methods. It provides the benefit of Godunov methods in capturing high Lorentz boosted flows but without complicated Riemann solvers, and the advantages of traditional artificial viscosity methods in their speed and flexibility. Additionally, the GRMHD equations are solved on an unstructured grid that supports local adaptive mesh refinement using a fully threaded oct-tree (in three dimensions) network to traverse the grid hierarchy across levels and immediate neighbors. A number of tests are presented to demonstrate robustness of the numerical algorithms and adaptive mesh framework over a wide spectrum of problems, boosts, and astrophysical applications, including relativistic shock tubes, shock collisions, magnetosonic shocks, Alfven wave propagation, blast waves, magnetized Bondi flow, and the magneto-rotational instability in Kerr black hole spacetimes.

A high-order accurate flow solution method for quadrilateral elements is developed and tested. The method employs the flux reconstruction (FR) approach proposed by Huynh. Ideal orders of accuracy are achieved for several problems on structured-based grids while some deterioration is observed on purely unstructured grids. Cares for curved wall boundary are also shown. The localized artificial diffusivity (LAD) method is incorporated to capture discontinuities that is crucial for compressible flow simulations. The FR+LAD works very well for 1D test cases but multi-dimensional shock problems requires further modifications to confine the artificial diffusivity to the shock location. Hybrid scheme of the FR and finite difference method taking advantage of 1D-feature alleviates the instability. I. Introduction Expanding demands for the computational fluid dynamics (CFD) requires more reliable numerical schemes. Specifically, turbulent flow simulations, acoustic simulation, combustion simulations and their combined problems. The key words in such simulations may be high accuracy and high resolution. In structured grid CFD, a compact finite difference scheme satisfies both requirements at higher level than existing other methods. If we put more emphasis on easier handling of complex geometries and flow-data dependent grid generation (grid adaptation), then unstructured grid method should also be explored and it should depart from 2nd-order scheme with limiter. Recent efforts in this area are development of the discontinuous Galerkin (DG) method, 1,2 the spectral finite volume (SV) method 3 and the staggered-grid multi-domain method spectral finite difference (SG) method 4,5 and its evolutional method called a spectral difference method. 6 The common and maybe the only way for those unstructured methods is to add degrees of freedom in each cell to achieve high order since the wide stencil of the original grid is not possible due to unstructured nature. The ENO-type method that uses original degrees of freedom in the grids and extends adaptive stencils to the flow solutions is also explored but the difficulty of extending stencils due to undirectional nature of unstructured grids and the unsuitability for parallel computing prohibits wide applications. A flux reconstruction (FR) scheme 7,8 studied in the paper is the new one among the former group. It uses differential form of conservation laws like the SD method but it is more general and includes equivalent counter part of the schemes in the group. In this paper, the performance of the FR is tested by a newly developed Navier-Stokes solver. Realistic flow applications of the FR in the Navier-Stokes equations are seldom seen since it is relatively new scheme. The Localized artificial diffusivity (LAD) 9–13 is incorporated for the shock capturing in the FR. The last item is motivated by the challenge and success of combining the SD and LAD in the recent study. 11,12 Wang 14 developed the lifting collocation penalty (LCP) method to extend the FR to triangular or tetrahedral grids and applied the LCP-FR combined method for hybrid unstructured grids. In the present study, only quadrilateral cells are treated to directly apply the LAD method that was developed on structured grids. The shock capturing by the FR+LAD high order scheme is tried.

A coupled regional and local model is required when groundwater flow and solute transport are to be simulated in local areas of interest with a finer grid while regional aquifer boundary and major stresses should be retained with a coarser grid. The coupled model should also maintain interactions between the regional and local flow systems. In the Beijing Plain (China), assessment of managed aquifer recharge (MAR), groundwater pollution caused by rivers, capture zone of well fields, and land subsidence at the cone of depression requires a coupled regional and local model. This study evaluates three methods for coupling regional and local flow models for simulating MAR in the Chaobai River catchment in the Beijing Plain. These methods are the conventional grid refinement (CGR) method, the local grid refinement (LGR) method and the unstructured grid (USG) method. The assessment included the comparison of the complexity of the coupled model construction, the goodness of fit of the computed and observed groundwater heads, the consistency of regional and local groundwater budgets, and the capture zone of a well filed influenced by the MAR site. The results indicated that the CGR method based on MODFLOW-2005 is the easiest to implement the coupled model, capable of reproducing regional and local groundwater heads and budget, and already coupled with density and viscosity dependent model codes for transport simulation. However, the CGR method inherits shortcomings of finite difference grids to create multiple local models with inefficient computing efforts. The USG method based on MODFLOW-USG has the advantage of creating multi-scale models and is flexible to simulate rivers, wells, irregular boundaries, heterogeneities and the MAR site. However, it is more difficult to construct the coupled models with the unstructured grids, therefore, a good graphic user interface is necessary for efficient model construction. The LGR method based on MODFLOW-LGR can be used to create multiple local models in uniform aquifer systems. So far, little effort has been devoted to upgrade the LGR method for complex aquifer structures and develop the coupled transport models.

A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids

The edge based Galerkin finite element formulation is used as the basic building block for the construction of multidimensional generalizations, on unstructured grids, of several higher order upwind biased procedures originally designed for the solution of the 1D compressible Euler system of equations. The use of a central type discretization for the viscous flux terms enables the simulation of multidimensional flows governed by the laminar compressible Navier Stokes equations. Numerical issues related to the development and implementation of multidimensional solution algorithms are considered. A number of inviscid and viscous flow simulations, in different flow regimes, are analyzed to enable the reader to assess the performance of the surveyed formulations.

In this paper, we focus on the transient convection problems of which physically bounds are already known. In tradition, the high‐resolutions (HR) schemes that use some parameters to measure the local monotonicity in the predicted scalar field, including the well‐known TVD‐type and NVD‐type schemes, are widely used. Recently, Qin Qian, et al., [1] proposed a new approach which use the upwind approximation as the ‘unbounded indicator”. During the computation, the actual scheme switches between QUICK scheme and 1st‐order upwind (FUD) scheme, according to whether the indicator locates in the range of physical bounds or not. However, the new approach does not apply to unstructured grids as the far‐upwind point for QUICK scheme is not easy to determine. Considering the easy implementation of central differencing (CD) scheme on structured and unstructured meshes, we propose a new convective scheme following the idea of [1], but using CD scheme as the high‐order base scheme. The new scheme is extremely simple, but with 2nd‐order accuracy. And best of all, the scheme can be applied to arbitrary grids. Testing is carried out on the 1D (transport of a top‐hat pulse) problem on a uniform grid. Further testing on the problem of advection of a square‐shaped scalar field on structured and unstructured grids are also done. Finally, we compare the solution obtained by FUD, Gamma, MUSCL and the proposed scheme and find that the new scheme performs well for all the test cases considered.

This paper is an overview of the finite element methods developed by the Team for Advanced Flow Simulation and Modeling (T*AFSM) [http://www.mems.rice.edu/TAFSM/] for computation of flow problems with moving boundaries and interfaces. This class of problems include those with free surfaces, two-fluid interfaces, fluid-object and fluid-structure interactions, and moving mechanical components. The methods developed can be classified into two main categories. The interface-tracking methods are based on the Deforming-Spatial-Domain/Stabilized Space-Time (DSD/SST) formulation, where the mesh moves to track the interface, with special attention paid to reducing the frequency of remeshing. The interface-capturing methods, typically used for free-surface and two-fluid flows, are based on the stabilized formulation, over non-moving meshes, of both the flow equations and the advection equation governing the time-evolution of an interface function marking the location of the interface. In this category, when it becomes neccessary to increase the accuracy in representing the interface beyond the accuracy provided by the existing mesh resolution around the interface, the Enhanced-Discretization Interface-Capturing Technique (EDICT) can be used to to accomplish that goal. In development of these two classes of methods, we had to keep in mind the requirement that the methods need to be applicable to 3D problems with complex geometries and that the associated large-scale computations need to be carried out on parallel computing platforms. Therefore our parallel implementations of these methods are based on unstructured grids and on both the distributed and shared memory parallel computing approaches. In addition to these two main classes of methods, a number of other ideas and methods have been developed to increase the scope and accuracy of these two classes of methods. The review of all these methods in our presentation here is supplemented by a number numerical examples from parallel computation of complex, 3D flow problems.

A parallel finite element method for incompressible Navier–Stokes flows based on unstructured grids

The application of unstructured grids can improve the solution of total field electromagnetic problems because these grids allow efficient local refinement of the mesh at the locations of high field gradients. Unstructured grids also provide the flexibility required for representing arbitrary topography and subsurface interfaces. We investigated the generalization of the standard Yee’s staggered scheme to unstructured tetrahedral-Voronoï grids using a finite-volume approach. We discretized the Helmholtz equation for the electric field in the frequency domain and solved the problem to find the projection of the total electric field along the edges of the tetrahedral elements. To compute the electric and magnetic fields at the observation points, an interpolation technique was employed, which uses the edge vector interpolation functions of the tetrahedral elements. To verify the presented scheme, two examples were evaluated, which revealed the computation of the total and secondary fields due to electric and magnetic sources in half-spaces that contain anomalous bodies. The results had good agreement with those from the literature. For the second example, accuracy studies were conducted to better understand the relative importance of different refinements in the grid and also the effect of the general quality of the grid. The results revealed the utmost importance of the quality of the tetrahedral grid and refinement at the observation locations. To validate the versatility of the approach, we synthesized helicopter-borne data for a model of the Ovoid ore body at Voisey’s Bay, Labrador, Canada, which had good agreement with real data.