Search

In the context of phase‐field modeling of fractures in incompressible materials, a mixed form of the elasticity equation can overcome possible volume locking effects. The drawback is that a coupled variational inequality system with three unknowns (displacements, pressure and phase‐field) has to be solved, which increases the overall workload. Efficient preconditioning at this point is an indispensable tool. In this work, a problem‐specific iterative solver is proposed leveraging the saddle‐point structure of the displacement and pressure variable. A Schur‐type preconditioner is developed to avoid ill‐conditioning of the phase‐field fracture problem. Finally, we show numerical results of a pressure‐driven benchmark which to confirm the robustness of the solver. [ABSTRACT FROM AUTHOR]

Problems arising in Earth's mantle convection involve finding the solution to Stokes systems with large viscosity contrasts. These systems contain localized features which, even with adaptive mesh refinement, result in linear systems that can be on the order of 109 or more unknowns. One common approach for preconditioning to the velocity block of these systems is to apply an Algebraic Multigrid (AMG) V‐cycle (as is done in the ASPECT software, for example), however, we find that AMG is lacking robustness with respect to problem size and number of parallel processes. Additionally, we see an increase in iteration counts with refinement when using AMG. In contrast, the Geometric Multigrid (GMG) method, by using information about the geometry of the problem, should offer a more robust option.Here we present a matrix‐free GMG V‐cycle which works on adaptively refined, distributed meshes, and we will compare it against the current AMG preconditioner (Trilinos ML) used in the ASPECT1 software. We will demonstrate the robustness of GMG with respect to problem size and show scaling up to 114,688 cores and 217 billion unknowns. All computations are run using the open‐source, finite element library deal.II.2 [ABSTRACT FROM AUTHOR]

We present the design and implementation details of a geometric multigrid method on adaptively refined meshes for massively parallel computations. The method uses local smoothing on the refined part of the mesh. Partitioning is achieved by using a space filling curve for the leaf mesh and distributing ancestors in the hierarchy based on the leaves. We present a model of the efficiency of mesh hierarchy distribution and compare its predictions to runtime measurements. The algorithm is implemented as part of the deal.II finite-element library and as such available to the public. [ABSTRACT FROM AUTHOR]

We study discretizations of the incompressible Navier–Stokes equations, written in the newly developed energy–momentum–angular momentum conserving (EMAC) formulation. We consider linearizations of the problem, which at each time step will reduce the computational cost, but can alter the conservation properties. We show that a skew-symmetrized linearization delivers the correct balance of (only) energy and that the Newton linearization conserves momentum and angular momentum, but conserves energy only up to the nonlinear residual. Numerical tests show that linearizing with 2 Newton steps at each time step is very effective at preserving all conservation laws at once, and giving accurate answers on long time intervals. The tests also show that the skew-symmetrized linearization is significantly less accurate. The tests also show that the Newton linearization of EMAC finite element formulation compares favorably to other traditionally used finite element formulation of the incompressible Navier–Stokes equations in primitive variables. [ABSTRACT FROM AUTHOR]

We present a scalable, parallel implementation of a solver for the solution of a phase‐field model for quasi‐static brittle fracture. The code is available as open source. Numerical solutions in 2d and 3d with adaptive mesh refinement show optimal scaling of the linear solver based on algebraic multigrid, and convergence of the phase‐field model towards exact values of functionals of interests such as the crack opening displacement or the total crack volume. In contrast to uniform refinement, adaptive mesh refinement allows us to recover optimal convergence rates for the non‐smooth solutions encountered in typical test problems. We also present numerical studies of the influence of the finite domain size on functional evaluations used to approximate the infinite domain. [ABSTRACT FROM AUTHOR]

In this paper, we propose a new numerical method to find the ground state energy of a given physical system within the Kohn–Sham density functional theory. The h‐adaptive finite element method is adopted for spatial discretization and implemented with matrix‐free operator evaluation. The ground state energy is found by performing unconstrained minimization with non‐orthogonal orbitals using the limited memory Broyden–Fletcher–Goldfarb–Shanno (BFGS) method. A geometric multigrid preconditioner is applied to improve the convergence. The clear advantage of the proposed approach is demonstrated on selected examples by comparing the performance to other methods such as preconditioned steepest descent minimization. The proposed method provides a solid framework toward O(N) complexity for the locally adaptive real‐space solution of density functional theory with finite elements. [ABSTRACT FROM AUTHOR]

Computations have helped elucidate the dynamics of Earth's mantle for several decades already. The numerical methods that underlie these simulations have greatly evolved within this time span, and today include dynamically changing and adaptively refined meshes, sophisticated and efficient solvers, and parallelization to large clusters of computers. At the same time, many of the methods - discussed in detail in a previous paper in this series - were developed and tested primarily using model problems that lack many of the complexities that are common to the realistic models our community wants to solve today. With several years of experience solving complex and realistic models, we here revisit some of the algorithm designs of the earlier paper and discuss the incorporation of more complex physics. In particular, we re-consider time stepping and mesh refinement algorithms, evaluate approaches to incorporate compressibility, and discuss dealing with strongly varying material coefficients, latent heat, and how to track chemical compositions and heterogeneities. Taken together and implemented in a high-performance, massively parallel code, the techniques discussed in this paper then allow for high resolution, 3-D, compressible, global mantle convection simulations with phase transitions, strongly temperature dependent viscosity and realistic material properties based on mineral physics data. [ABSTRACT FROM AUTHOR]

We study conservation properties of Galerkin methods for the incompressible Navier–Stokes equations, without the divergence constraint strongly enforced. In typical discretizations such as the mixed finite element method, the conservation of mass is enforced only weakly, and this leads to discrete solutions which may not conserve energy, momentum, angular momentum, helicity, or vorticity, even though the physics of the Navier–Stokes equations dictate that they should. We aim in this work to construct discrete formulations that conserve as many physical laws as possible without utilizing a strong enforcement of the divergence constraint, and doing so leads us to a new formulation that conserves each of energy, momentum, angular momentum, enstrophy in 2D, helicity and vorticity (for reference, the usual convective formulation does not conserve most of these quantities). Several numerical experiments are performed, which verify the theory and test the new formulation. [ABSTRACT FROM AUTHOR]

Melt generation and migration are an important link between surface processes and the thermal and chemical evolution of the Earth's interior. However, their vastly different timescales make it difficult to study mantle convection and melt migration in a unified framework, especially for 3-D global models. And although experiments suggest an increase in melt volume of up to 20 per cent from the depth of melt generation to the surface, previous computations have neglected the individual compressibilities of the solid and the fluid phase. Here, we describe our extension of the finite element mantle convection code ASPECT that adds melt generation and migration. We use the original compressible formulation of the McKenzie equations, augmented by an equation for the conservation of energy. Applying adaptive mesh refinement to this type of problems is particularly advantageous, as the resolution can be increased in areas where melt is present and viscosity gradients are high, whereas a lower resolution is sufficient in regions without melt. Together with a high-performance, massively parallel implementation, this allows for high-resolution, 3-D, compressible, global mantle convection simulations coupled with melt migration. We evaluate the functionality and potential of this method using a series of benchmarks and model setups, compare results of the compressible and incompressible formulation, and show the effectiveness of adaptive mesh refinement when applied to melt migration. Our model of magma dynamics provides a framework for modelling processes on different scales and investigating links between processes occurring in the deep mantle and melt generation and migration. This approach could prove particularly useful applied to modelling the generation of komatiites or other melts originating in greater depths. The implementation is available in the Open Source ASPECT repository. [ABSTRACT FROM AUTHOR]

We investigate a flux-preserving enforcement of inhomogeneous Dirichlet boundary conditions for velocity, u | ∂ Ω = g , for use with finite element methods for incompressible flow problems that strongly enforce mass conservation. Typical enforcement via nodal interpolation is not flux-preserving in general, and it can create divergence error even when divergence-free elements are used. We show with analysis and numerical tests that by slightly (and locally) changing nodal interpolation, the enforcement recovers flux-preservation, is optimally accurate, and delivers divergence-free solutions when used with divergence-free finite elements. [ABSTRACT FROM AUTHOR]

We derive a new kind of boundary conditions for the vorticity equation with solid wall boundaries for fluid flow problems. The formulation uses a Dirichlet condition for the normal component of vorticity and Neumann type conditions for the tangential components. In a Galerkin (integral) formulation the tangential condition is natural, i.e., it is enforced by a right-hand side functional and does not impose a boundary constraint on trial and test spaces. The functional involves the pressure variable, and we discuss several velocity–vorticity formulations where the proposed condition is appropriate. Several numerical experiments are given that illustrate the validity of the new approach. [ABSTRACT FROM AUTHOR]

In this paper, we consider phase-field based fracture propagation in elastic media. The main purpose is the development of a robust and efficient numerical scheme. To enforce crack irreversibility as a constraint, we use a primal-dual active set strategy, which can be identified as a semi-smooth Newton method. The active set iteration is merged with the Newton iteration for solving the fully-coupled nonlinear partial differential equation discretized using finite elements, resulting in a single, rapidly converging nonlinear scheme. It is well known that phase-field models require fine meshes to accurately capture the propagation dynamics of the crack. Because traditional estimators based on adaptive mesh refinement schemes are not appropriate, we develop a predictor-corrector scheme for local mesh adaptivity to reduce the computational cost. This method is both robust and efficient and allows us to treat temporal and spatial refinements and to study the influence of model regularization parameters. Finally, our proposed approach is substantiated with different numerical tests for crack propagation in elastic media and pressurized fracture propagation in homogeneous and heterogeneous media. [ABSTRACT FROM AUTHOR]

SUMMARY Numerical simulation of the processes in the Earth's mantle is a key piece in understanding its dynamics, composition, history and interaction with the lithosphere and the Earth's core. However, doing so presents many practical difficulties related to the numerical methods that can accurately represent these processes at relevant scales. This paper presents an overview of the state of the art in algorithms for high-Rayleigh number flows such as those in the Earth's mantle, and discusses their implementation in the Open Source code Aspect (Advanced Solver for Problems in Earth's ConvecTion). Specifically, we show how an interconnected set of methods for adaptive mesh refinement (AMR), higher order spatial and temporal discretizations, advection stabilization and efficient linear solvers can provide high accuracy at a numerical cost unachievable with traditional methods, and how these methods can be designed in a way so that they scale to large numbers of processors on compute clusters. Aspect relies on the numerical software packages deal.II and Trilinos, enabling us to focus on high level code and keeping our implementation compact. We present results from validation tests using widely used benchmarks for our code, as well as scaling results from parallel runs. [ABSTRACT FROM AUTHOR]

This paper provides an overview of the new features of the finite element library deal.II, version 9.5. [ABSTRACT FROM AUTHOR]

Mantle convection models based on geophysical constraints have provided us with a basic understanding of the forces driving and resisting plate motions on Earth. However, existing studies computing the balance of underlying forces are contradicting, and the impact of plate boundary geometry on surface deformation remains unknown. We address these issues by developing global instantaneous 3‐D mantle convection models with a heterogeneous density and viscosity distribution and weak plate boundaries prescribed using different geometries. We find that the plate boundary geometry of the Global Earthquake Model (GEM, Pagani et al., 2018, https://doi.org/10.1177/8755293020931866), featuring open plate boundaries with discrete lithospheric‐depth weak zones in the oceans and distributed crustal faults within continents, achieves the best fit to the observed GPS data with a directional correlation of 95.1% and a global point‐wise velocity residual of 1.87 cm/year. A good fit also requires plate boundaries being 3 to 4 orders of magnitude weaker than the surrounding lithosphere and low asthenospheric viscosities between 5 × 1017 and 5 × 1018 Pa s. Models without asthenospheric and lower mantle heterogeneities retain on average 30% and 70% of the plate speeds, respectively. Our results show that Earth's plate boundaries are not uniform and better described by more discrete plate boundaries within the oceans and distributed faults within continents. Furthermore, they emphasize the impact of plate boundary geometry on the direction and speed of plate motions and reaffirm the importance of slab pull in the uppermost mantle as a major plate driving force. Plain Language Summary: Plate tectonics can explain several geological and geophysical phenomena on Earth and is closely coupled to convection in the underlying mantle. To understand this plate–mantle coupling and quantify the forces contributing to plate motion, we develop high‐resolution three‐dimensional computational models of the Earth's present‐day mantle flow utilizing available geophysical constraints on density distribution and rheology. Additionally, we prescribe weak zones at the location of plate boundaries. We use different plate boundary geometries, forming either open or closed polygons, and we vary how easily the plate boundaries and the asthenosphere directly below the plates can be deformed to determine which model best fits observed plate motions. Our best‐fitting model features open plate boundaries that are weak (∼4 orders of magnitude weaker than the surrounding lithosphere) and traverse the whole plate in the oceans, but are shallower and more distributed within continents. The asthenosphere in these models is even weaker than the plate boundaries. Furthermore, we find that the downward force caused by subducted slabs contributes the most to the observed surface velocities. Our models suggest that plate boundaries are not uniformly weak everywhere and that their geometry has a strong influence on the direction and speed of plate motion. Key Points: We model plate motions in global instantaneous 3‐D mantle convection models with different plate boundary geometriesEarth's plate boundaries are not uniform and better described by discrete shear zones in the oceans and distributed faults within continentsSlab pull within the uppermost mantle (<300 km depth) contributes about 70% of the total plate driving force [ABSTRACT FROM AUTHOR]

Due to the increasing availability of high-performance computing over the past few decades, numerical models have become an important tool for research in geodynamics. Several generations of mantle convection software have been developed, but due to their differing methods and increasing complexity it is important to evaluate the accuracy of each new model generation to ensure published geodynamic research is reliable and reproducible. Here we explore the accuracy of the open-source, finite-element codes ASPECT and CitcomS as a function of mesh spacing using low to moderate-Rayleigh-number models in steady-state thermal convection. ASPECT (Advanced Solver for Problems in Earth's ConvecTion) is a new-generation mantle convection code that enables modeling global mantle convection with realistic parameters and complicated physical processes using adaptive mesh refinement. We compare the ASPECT results with calculations from the finite-element code CitcomS , which has a long history of use in the geodynamics community. We find that the globally averaged quantities, i.e., root-mean-square (rms) velocity, mean temperature, and Nusselt number at the top and bottom of the shell, agree to within 1 % (and often much better) for calculations with sufficient mesh resolution. We also show that there is excellent agreement of the time evolution of both the rms velocity and the Nusselt numbers between the two codes for otherwise identical parameters. Based on our results, we are optimistic that similar agreement would be achieved for calculations performed at the convective vigor expected for Earth, Venus, and Mars. [ABSTRACT FROM AUTHOR]

Geodynamic simulations increasingly rely on models with a true free surface to investigate questions of dynamic topography, tectonic deformation, gravity perturbations, and global mantle convection. However, implementations of free surface boundary conditions have proven challenging from a standpoint of accuracy, robustness, and stability. In particular, time integration of a free surface tends to suffer from a numerical instability that manifests as sloshing surface motions, also known as the “drunken sailor” instability. This instability severely limits stable timestep sizes to those much smaller than can be used in geodynamic simulations without a free surface. Several schemes have been proposed in the literature to deal with these instabilities. Here we analyze the problem of creeping viscous flow with a free surface and discuss the origin of these instabilities. We demonstrate their cause and how existing stabilization schemes work to damp them out. We also propose a new scheme for removing instabilities from free surface calculations. It does not require modifications to the system matrix, nor additional variables, but is instead an explicit scheme based on nonstandard finite differences. It relies on a single stabilization parameter which may be identified with the smallest relaxation timescale of the free surface. Finally, we present numerical results to show the effectiveness of the new approach and discuss the free surface implementation in the open source, community based mantle convection software ASPECT . [ABSTRACT FROM AUTHOR]

This paper provides an overview of the new features of the finite element library deal.II version 8.4. [ABSTRACT FROM AUTHOR]

Due to the increasing availability of high-performance computing over the past decades numerical models have become an important tool for research in geodynamics. Several generations of mantle convection software have been developed, but due to their differing methods and increasing complexity it is important to evaluate the accuracy of each new model generation to ensure published geodynamic research is reliable and reproducible. We here explore the accuracy of the open-source, finite-element codes ASPECT and CitcomS as a function of mesh spacing using low to moderate Rayleigh number models in steady-state, thermal convection. ASPECT (Advanced Solver for Problems in Earth's ConvecTion) is a new generation mantle convection code that enables modeling global mantle convection with realistic parameters and complicated physical processes using adaptive mesh refinement (Kronbichler et al., 2012; Heister et al., 2017). We compare the ASPECT results with calculations from the finite element code CitcomS (Zhong et al., 2000; Tan et al., 2006; Zhong et al., 2008), which has a long history of use in the geodynamics community. We find that the globally-averaged quantities: RMS velocity, mean temperature, and Nusselt number at the top and bottom of the shell, agree to within 1%, and often much better, for calculations with sufficient mesh resolution. We also show that there is excellent agreement of the time-evolution of both the RMS velocity and the Nusselt numbers between the two codes for otherwise identical parameters. Based on our results we are optimistic that similar agreement would be achieved for calculations performed at convective vigor expected for Earth, Venus, and Mars. [ABSTRACT FROM AUTHOR]

This paper provides an overview of the new features of the finite element library deal.II, version 9.4. [ABSTRACT FROM AUTHOR]

Abstract: In this paper the grad–div stabilization for the incompressible Navier–Stokes finite element approximations is considered from two different viewpoints: (i) as a variational multiscale approach for the pressure subgrid modeling and (ii) as a stabilization procedure of least-square type. Some new error estimates for the linearized problem with the grad–div stabilization are proved with the help of norms induced by the pressure Schur complement operator. We discuss the stabilization parameter choice arising in the frameworks of least-square and multiscale methods and consider assumptions which allow to relate both approaches. [Copyright &y& Elsevier]

This paper provides an overview of the new features of the finite element library deal.II, version 9.3. [ABSTRACT FROM AUTHOR]

deal.II is a state-of-the-art finite element library focused on generality, dimension-independent programming, parallelism, and extensibility. Herein, we outline its primary design considerations and its sophisticated features such as distributed meshes, h p -adaptivity, support for complex geometries, and matrix-free algorithms. But deal.II is more than just a software library: It is also a diverse and worldwide community of developers and users, as well as an educational platform. We therefore also discuss some of the technical and social challenges and lessons learned in running a large community software project over the course of two decades. [ABSTRACT FROM AUTHOR]

This paper provides an overview of the new features of the finite element library deal.II, version 9.2. [ABSTRACT FROM AUTHOR]

Mantle convection and long-term lithosphere dynamics in the Earth and other planets can be treated as the slow deformation of a highly viscous fluid, and as such can be described using the compressible Navier–Stokes equations. Since on Earth-sized planets the influence of compressibility is not a dominant effect, density deviations from a reference profile are at most on the order of a few percent and using the full governing equations poses numerical challenges, most modelling studies have simplified the governing equations. Common approximations assume a temporally constant, but depth-dependent reference profile for the density (the anelastic liquid approximation), or drop compressibility altogether and use a constant reference density (the Boussinesq approximation). In most previous studies of mantle convection and crustal dynamics, one can assume that the error introduced by these approximations was small compared to the errors that resulted from poorly constrained material behaviour and limited numerical accuracy. However, as model parametrizations have become more realistic, and model resolution has improved, this may no longer be the case and the error due to using simplified conservation equations might no longer be negligible: while such approximations may be reasonable for models of mantle plumes or slabs traversing the whole mantle, they may be unsatisfactory for layered materials experiencing phase transitions or materials undergoing significant heating or cooling. For example, at boundary layers or close to dynamically changing density gradients, the error arising from the use of the aforementioned compressibility approximations can be the dominant error source, and common approximations may fail to capture the physical behaviour of interest. In this paper, we discuss new formulations of the continuity equation that include dynamic density variations due to temperature, pressure and composition without using a reference profile for the density. We quantify the improvement in accuracy relative to existing formulations in a number of benchmark models and evaluate for which practical applications these effects are important. Finally, we consider numerical aspects of the new formulations. We implement and test these formulations in the freely available community software aspect , and use this code for our numerical experiments. [ABSTRACT FROM AUTHOR]

This paper provides an overview of the new features of the finite element library deal.II, version 9.1. [ABSTRACT FROM AUTHOR]

Modelling the porous flow of melt through a viscously deforming solid rock matrix is a useful tool for interpreting observations from the Earth's surface, and advances our understanding of the dynamics of the Earth's interior. However, the system of equations describing this process becomes mathematically degenerate in the limit of vanishing melt fraction. Numerical methods that do not consider this degeneracy or avoid it solely by regularizing specific material properties generally become computationally expensive as soon as the melt fraction approaches zero in some part of the domain. Here, we present a new formulation of the equations for coupled magma/mantle dynamics that addresses this problem, and allows it to accurately compute large-scale 3-D magma/mantle dynamics simulations with extensive regions of zero melt fraction. We achieve this by rescaling one of the solution variables, the compaction pressure, which ensures that for vanishing melt fraction, the equation causing the degeneracy becomes an identity and the other two equations revert to the Stokes system. This allows us to split the domain into two parts: in mesh cells where melt is present, we solve the coupled system of magma/mantle dynamics. In cells without melt, we solve the Stokes system as it is done for mantle convection without melt transport and constrain the remaining degrees of freedom. We have implemented this formulation in the open source geodynamic modelling code Aspect and illustrate the improved performance compared to the previous three-field formulation, showing numerically that the new formulation is robust in terms of problem size and only slightly sensitive to model parameters. Beyond that, we demonstrate the applicability to realistic problems by showing large-scale 2-D and 3-D models of mid-ocean ridges with complex rheology. Hence, we believe that our new formulation and its implementation in Aspect will prove a valuable tool for studying the interaction of melt segregating through and interacting with a solid host rock in the Earth and other planetary bodies using high-resolution, 3-D simulations. [ABSTRACT FROM AUTHOR]

This paper provides an overview of the new features of the finite element library deal.II version 9.0. [ABSTRACT FROM AUTHOR]

This paper provides an overview of the new features of the finite element library deal.II version 8.5. [ABSTRACT FROM AUTHOR]