skip to main content


Title: Evaluating the accuracy of hybrid finite element/particle-in-cell methods for modelling incompressible Stokes flow
SUMMARY

Combining finite element methods for the incompressible Stokes equations with particle-in-cell methods is an important technique in computational geodynamics that has been widely applied in mantle convection, lithosphere dynamics and crustal-scale modelling. In these applications, particles are used to transport along properties of the medium such as the temperature, chemical compositions or other material properties; the particle methods are therefore used to reduce the advection equation to an ordinary differential equation for each particle, resulting in a problem that is simpler to solve than the original equation for which stabilization techniques are necessary to avoid oscillations.

On the other hand, replacing field-based descriptions by quantities only defined at the locations of particles introduces numerical errors. These errors have previously been investigated, but a complete understanding from both the theoretical and practical sides was so far lacking. In addition, we are not aware of systematic guidance regarding the question of how many particles one needs to choose per mesh cell to achieve a certain accuracy.

In this paper we modify two existing instantaneous benchmarks and present two new analytic benchmarks for time-dependent incompressible Stokes flow in order to compare the convergence rate and accuracy of various combinations of finite elements, particle advection and particle interpolation methods. Using these benchmarks, we find that in order to retain the optimal accuracy of the finite element formulation, one needs to use a sufficiently accurate particle interpolation algorithm. Additionally, we observe and explain that for our higher-order finite-element methods it is necessary to increase the number of particles per cell as the mesh resolution increases (i.e. as the grid cell size decreases) to avoid a reduction in convergence order.

Our methods and results allow designing new particle-in-cell methods with specific convergence rates, and also provide guidance for the choice of common building blocks and parameters such as the number of particles per cell. In addition, our new time-dependent benchmark provides a simple test that can be used to compare different implementations, algorithms and for the assessment of new numerical methods for particle interpolation and advection. We provide a reference implementation of this benchmark in aspect (the ‘Advanced Solver for Problems in Earth’s ConvecTion’), an open source code for geodynamic modelling.

 
more » « less
Award ID(s):
1835673
NSF-PAR ID:
10120480
Author(s) / Creator(s):
 ;  ;  ;  
Publisher / Repository:
Oxford University Press
Date Published:
Journal Name:
Geophysical Journal International
Volume:
219
Issue:
3
ISSN:
0956-540X
Page Range / eLocation ID:
p. 1915-1938
Format(s):
Medium: X
Sponsoring Org:
National Science Foundation
More Like this
  1. Simulation of flow and transport in petroleum reservoirs involves solving coupled systems of advection-diffusion-reaction equations with nonlinear flux functions, diffusion coefficients, and reactions/wells. It is important to develop numerical schemes that can approximate all three processes at once, and to high order, so that the physics can be well resolved. In this paper, we propose an approach based on high order, finite volume, implicit, Weighted Essentially NonOscillatory (iWENO) schemes. The resulting schemes are locally mass conservative and, being implicit, suited to systems of advection-diffusion-reaction equations. Moreover, our approach gives unconditionally L-stable schemes for smooth solutions to the linear advection-diffusion-reaction equation in the sense of a von Neumann stability analysis. To illustrate our approach, we develop a third order iWENO scheme for the saturation equation of two-phase flow in porous media in two space dimensions. The keys to high order accuracy are to use WENO reconstruction in space (which handles shocks and steep fronts) combined with a two-stage Radau-IIA Runge-Kutta time integrator. The saturation is approximated by its averages over the mesh elements at the current time level and at two future time levels; therefore, the scheme uses two unknowns per grid block per variable, independent of the spatial dimension. This makes the scheme fairly computationally efficient, both because reconstructions make use of local information that can fit in cache memory, and because the global system has about as small a number of degrees of freedom as possible. The scheme is relatively simple to implement, high order accurate, maintains local mass conservation, applies to general computational meshes, and appears to be robust. Preliminary computational tests show the potential of the scheme to handle advection-diffusion-reaction processes on meshes of quadrilateral gridblocks, and to do so to high order accuracy using relatively long time steps. The new scheme can be viewed as a generalization of standard cell-centered finite volume (or finite difference) methods. It achieves high order in both space and time, and it incorporates WENO slope limiting. 
    more » « less
  2. Abstract

    We develop a generalized interpolation material point method (GIMPM) for the shallow shelf approximation (SSA) of ice flow. The GIMPM, which can be viewed as a particle version of the finite element method, is used here to solve the shallow shelf approximations of the momentum balance and ice thickness evolution equations. We introduce novel numerical schemes for particle splitting and integration at domain boundaries to accurately simulate the spreading of an ice shelf. The advantages of the proposed GIMPM‐SSA framework include efficient advection of history or internal state variables without diffusion errors, automated tracking of the ice front and grounding line at sub‐element scales, and a weak formulation based on well‐established conventions of the finite element method with minimal additional computational cost. We demonstrate the numerical accuracy and stability of the GIMPM using 1‐D and 2‐D benchmark examples. We also compare the accuracy of the GIMPM with the standard material point method (sMPM) and a reweighted form of the sMPM. We find that the grid‐crossing error is very severe for SSA simulations with the sMPM, whereas the GIMPM successfully mitigates this error. While the grid‐crossing error can be reasonably reduced in the sMPM by implementing a simple material point reweighting scheme, this approach it not as accurate as the GIMPM. Thus, we illustrate that the GIMPM‐SSA framework is viable for the simulation of ice sheet‐shelf evolution and enables boundary tracking and error‐free advection of history or state variables, such as ice thickness or damage.

     
    more » « less
  3. This paper develops a tree-topological local mesh refinement (TLMR) method on Cartesian grids for the simulation of bio-inspired flow with multiple moving objects. The TLMR nests refinement mesh blocks of structured grids to the target regions and arrange the blocks in a tree topology. The method solves the time-dependent incompressible flow using a fractional-step method and discretizes the Navier-Stokes equation using a finite-difference formulation with an immersed boundary method to resolve the complex boundaries. When iteratively solving the discretized equations across the coarse and fine TLMR blocks, for better accuracy and faster convergence, the momentum equation is solved on all blocks simultaneously, while the Poisson equation is solved recursively from the coarsest block to the finest ones. When the refined blocks of the same block are connected, the parallel Schwarz method is used to iteratively solve both the momentum and Poisson equations. Convergence studies show that the algorithm is second-order accurate in space for both velocity and pressure, and the developed mesh refinement technique is benchmarked and demonstrated by several canonical flow problems. The TLMR enables a fast solution to an incompressible flow problem with complex boundaries or multiple moving objects. Various bio-inspired flows of multiple moving objects show that the solver can save over 80% computational time, proportional to the grid reduction when refinement is applied. 
    more » « less
  4. Abstract

    In prior work we found that precise approximation of the continuity constraint is crucial for accurate propagation of tracer data when advected through a background incompressible velocity field (Sime et al., 2021,https://doi.org/10.1029/2020gc009349). Here we extend this investigation to compressible flows using the anelastic liquid approximation (ALA) and address four related issues: (a) Exact conservation of tracer discretized fields through a background compressible velocity; (b) Exact mass conservation; (c) Addition and removal of tracers without affecting (exact) conservation to preserve a consistent number of tracers per cell; and (d) the diffusion of tracer data, for example, as induced by thermal or chemical effects. In this process we also present an abstract formulation of the interior penalty hybrid discontinuous Galerkin (HDG) finite element formulation for diffusion problems and apply it to the advection‐diffusion and compressible Stokes systems. Finally we present numerical experiments exhibiting the HDG compressible Stokes momentum formulation's superconvergent compressibility approximation and reproduce examples of a community benchmark for the ALA.

     
    more » « less
  5. Abstract

    We introduce, analyze, and test an interpolation operator designed for use with continuous data assimilation (DA) of evolution equations that are discretized spatially with the finite element method. The interpolant is constructed as an approximation of theL2projection operator onto piecewise constant functions on a coarse mesh, but which allows nudging to be done completely at the linear algebraic level, independent of the rest of the discretization, with a diagonal matrix that is simple to construct; it can even completely remove the need for explicit construction of a coarse mesh. We prove the interpolation operator has sufficient stability and accuracy properties, and we apply it to algorithms for both fluid transport DA and incompressible Navier–Stokes DA. For both applications we prove the DA solutions with arbitrary initial conditions converge to the true solution (up to discretization error) exponentially fast in time, and are thus long‐time accurate. Results of several numerical tests are given, which both illustrate the theory and demonstrate its usefulness on practical problems.

     
    more » « less