Marcin Dabrowski

MILAMIN: MATLAB-based finite element method solver for large problems

The finite element method (FEM) combined with unstructured meshes forms an elegant and versatile approach capable of dealing with the complexities of problems in Earth science. Practical applications often require high-resolution models that necessitate advanced computational strategies. We therefore developed “Million a Minute” (MILAMIN), an efficient MATLAB implementation of FEM that is capable of setting up, solving, and postprocessing two-dimensional problems with one million unknowns in one minute on a modern desktop computer. MILAMIN allows the user to achieve numerical resolutions that are necessary to resolve the heterogeneous nature of geological materials. In this paper we provide the technical knowledge required to develop such models without the need to buy a commercial FEM package, programming compiler-language code, or hiring a computer specialist. It has been our special aim that all the components of MILAMIN perform efficiently, individually and as a package. While some of the components rely on readily available routines, we develop others from scratch and make sure that all of them work together efficiently. One of the main technical focuses of this paper is the optimization of the global matrix computations. The performance bottlenecks of the standard FEM algorithm are analyzed. An alternative approach is developed that sustains high performance for any system size. Applied optimizations eliminate Basic Linear Algebra Subprograms (BLAS) drawbacks when multiplying small matrices, reduce operation count and memory requirements when dealing with symmetric matrices, and increase data transfer efficiency by maximizing cache reuse. Applying loop interchange allows us to use BLAS on large matrices. In order to avoid unnecessary data transfers between RAM and CPU cache we introduce loop blocking. The optimization techniques are useful in many areas as demonstrated with our MILAMIN applications for thermal and incompressible flow (Stokes) problems. We use these to provide performance comparisons to other open source as well as commercial packages and find that MILAMIN is among the best performing solutions, in terms of both speed and memory usage. The corresponding MATLAB source code for the entire MILAMIN, including input generation, FEM solver, and postprocessing, is available from the authors (http://www.milamin.org) and can be downloaded as auxiliary material.

Parallel symmetric sparse matrix-vector product on scalar multi-core CPUs

We present a massively parallel implementation of symmetric sparse matrix-vector product for modern clusters with scalar multi-core CPUs. Matrices with highly variable structure and density arising from unstructured three-dimensional FEM discretizations of mechanical and diffusion problems are studied. A metric of the effective memory bandwidth is introduced to analyze the impact on performance of a set of simple, well-known optimizations: matrix reordering, manual prefetching, and blocking. A modification to the CRS storage improving the performance on multi-core Opterons is shown. The performance of an entire SMP blade rather than the per-core performance is optimized. Even for the simplest 4 node mechanical element our code utilizes close to 100% of the per-blade available memory bandwidth. We show that reducing the storage requirements for symmetric matrices results in roughly two times speedup. Blocking brings further storage savings and a proportional performance increase. Our results are compared to existing state-of-the-art implementations of SpMV, and to the dense BLAS2 performance. Parallel efficiency on 5400 Opteron cores of the Cray XT4 cluster is around 80-90% for problems with approximately 25^3 mesh nodes per core. For a problem with 820 million degrees of freedom the code runs with a sustained performance of 5.2 TeraFLOPs, over 20% of the theoretical peak.

Survival of LLSVPs for billions of years in a vigorously convecting mantle: Replenishment and destruction of chemical anomaly

We study segregation of the subducted oceanic crust (OC) at the core-mantle boundary and its ability to accumulate and form large thermochemical piles (such as the seismically observed Large Low Shear Velocity Provinces (LLSVPs)). Our high-resolution numerical simulations of thermochemical mantle convection suggest that the longevity of LLSVPs for up to three billion years, and possibly longer, can be ensured by a balance in the rate of segregation of high-density OC material to the core-mantle boundary (CMB) and the rate of its entrainment away from the CMB by mantle upwellings. For a range of parameters tested in this study, a large-scale compositional anomaly forms at the CMB, similar in shape and size to the LLSVPs. Neutrally buoyant thermochemical piles formed by mechanical stirring—where thermally induced negative density anomaly is balanced by the presence of a fraction of dense anomalous material—best resemble the geometry of LLSVPs. Such neutrally buoyant piles tend to emerge and survive for at least 3 Gyr in simulations with quite different parameters. We conclude that for a plausible range of values of density anomaly of OC material in the lower mantle—it is likely that it segregates to the CMB, gets mechanically mixed with the ambient material, and forms neutrally buoyant large-scale compositional anomalies similar in shape to the LLSVPs.

Monoclinic and triclinic 3D flanking structures around elliptical cracks

We use the Eshelby solution modified for a viscous fluid to model the evolution of three-dimensional flanking structures in monoclinic shear zones. Shearing of an elliptical crack strongly elongated perpendicular to the flow direction produces a cylindrical flanking structure which is reproducible with 2D plane strain models. In contrast, a circular or even narrow, slit-shaped crack exhibits a reduced magnitude of the velocity jump across the crack and results in smaller offset and a narrower zone of deflection than predicted with 2D-models. Even more significant deviations are observed if the crack axes are oriented at an oblique angle to the principal flow directions, where the velocity jump is oblique to the resolved shear direction and is modified during progressive deformation. The resulting triclinic geometry represents a rare example of triclinic structures developing in monoclinic flow and may be used to estimate the flow kinematics of the shear zone.

Efficient 3D stencil computations using CUDA

We present an efficient implementation of 7-point and 27-point stencils on high-end Nvidia GPUs. A new method of reading the data from the global memory to the shared memory of thread blocks is developed. The method avoids conditional statements and requires only two coalesced instructions to load the tile data with the halo (ghost zone). Additional optimizations include storing only one XY tile of data at a time in the shared memory to lower shared memory requirements, common subexpression elimination to reduce the number of instructions, and software prefetching to overlap arithmetic and memory instructions, and enhance latency hiding. The efficiency of our implementation is analyzed using a simple stencil memory footprint model that takes into account the actual halo overhead due to the minimum memory transaction size on the GPUs. Through experiments we demonstrate that in our implementation the memory overhead due to the halos is largely eliminated by good reuse of the halo data in the memory caches, and that our method of reading the data is close to optimal in terms of memory bandwidth usage. Detailed performance analysis for single precision stencil computations, and performance results for single and double precision arithmetic on two Tesla cards are presented. Our stencil implementations are more efficient than any other implementation described in the literature to date. On Tesla C2050 with single and double precision arithmetic our 7-point stencil achieves an average throughput of 12.3 and 6.5Gpts/s, respectively (98 GFLOP/s and 52 GFLOP/s, respectively). The symmetric 27-point stencil sustains a throughput of 10.9 and 5.8 Gpts/s, respectively.

An analytical benchmark with combined pressure and shear loading for elastoplastic numerical models

We discuss the benchmark strategy to check the accuracy of elastoplastic numerical solutions based on a fully two-dimensional analytical solution. Associated rate-independent non-hardening plasticity with von Mises or Tresca criteria is assumed throughout. Algorithms for integration of the rate equations, strategies for stress updating over a time step, and tangent operators are discussed. The accuracy of a simple incremental algorithm as a function of the time step is discussed.

Transient cluster formation in sheared non-Brownian suspensions.

We perform numerical simulations of non-Brownian suspensions in the laminar flow regime to study the scaling behavior of particle clusters and collisions under shear. As the particle fraction approaches the maximum packing fraction, large transient clusters appear in the system. We use methods from percolation theory to discuss the cluster size distribution. We also give a scaling relation for the percolation threshold as well as system size effects through time-dependent fluctuations of this threshold and relate them to system size. System size effects are important close to the maximum packing fraction due to the divergence of the cluster length scale. We then investigate the transient nature of the clusters through characterization of particle collisions and show that collision times exhibit scale-invariant properties. Finally, we show that particle collision times can be modeled as first-passage processes.

Nonlinear Viscoelastic Closure of Salt Cavities

Time-dependent hole closure is a major problem for the many cavities present in rock salt. We use analytical and numerical methods to study how cylindrical holes close under pressure loads with time. We treat salt as a viscoelastic fluid and we use an incompressible nonlinear Maxwell constitutive law to model its mechanical behavior. The viscosity is described by either a power law or an Ellis model depending on whether dislocation creep is considered alone or in combination with pressure solution. The instantaneous closure rate of a circular hole in a power law-based viscoelastic salt is fully determined analytically. A proxy for the transient closure velocity at the rim is also proposed based on a modified version of the characteristic relaxation time θ proposed by Wang et al. (Rock Mech Rock Eng 48(6):2369–2382, 2015) and it has less than 3% inaccuracy for times smaller than 3θ, irrespective of the load or salt type. We derive an analytical expression describing the instantaneous closure rate in an Ellis-based viscoelastic salt. A load threshold determines whether steady state is approached initially. The time θ is also a characteristic relaxation time for this constitutive law, and a master curve can be used to describe the evolution of the closure velocity with time. Using these characteristic values in a typical application underlines the importance of considering pressure solution, in addition to dislocation creep, when studying hole closure in rock salt.

Depth-averaged Lattice Boltzmann and Finite Element methods for single-phase flows in fractures with obstacles

Abstract We use Lattice Boltzmann Method (LBM) MRT and Cumulant schemes to study the performance and accuracy of single-phase flow modeling for propped fractures. The simulations are run using both the two- and three-dimensional Stokes equations, and a 2.5D Stokes–Brinkman approximate model. The LBM results are validated against Finite Element Method (FEM) simulations and an analytical solution to the Stokes–Brinkman flow around an isolated circular obstacle. Both LBM and FEM 2.5D Stokes–Brinkman models are able to reproduce the analytical solution for an isolated circular obstacle. In the case of 2D Stokes and 2.5D Stokes–Brinkman models, the differences between the extrapolated fracture permeabilities obtained with LBM and FEM simulations for fractures with multiple obstacles are below 1%. The differences between the fracture permeabilities computed using 3D Stokes LBM and FEM simulations are below 8% . The differences between the 3D Stokes and 2.5 Stokes–Brinkman results are less than 7% for FEM study, and 8% for the LBM case. The velocity perturbations that are introduced around the obstacles are not fully captured by the parabolic velocity profile inherent to the 2.5D Stokes–Brinkman model.

3D Fold Pattern Formation

Folds on all scales from millimeters to kilometers can be the result of the mechanical instability that arises when a mechanically stratified system is subjected to layer-parallel compression. While the resulting fold patterns are three dimensional, their geometries are often simplified by assuming that there is no shape variation in the third dimension. This facilitates the analysis and has resulted in a large number of studies that investigate the folding instability for a variety of rheologies: viscous, power-law, and anisotropic. Studies of three dimensional folding have mostly focused on analog models and geometric models. The latter led to the development of fold shape classification tools; in particular fold interference pattern classification (e.g., Grasemann et al. 2004, Odriscoll 1962, Ramsay 1967, Thiessen and Means 1980) and geometrical analysis based on differential geometry geometry (e.g., Lisle and Toimil 2007, Mynatt et al. 2007). The theoretical aspects of the mechanics of folding in three dimensional geological systems are only analyzed by in a few papers; Fletcher (1991, 1995), Ghosh (1970), Kaus and Schmalholz (2006), Muhlhaus (1998), and Schmid et al. (2008). The target of this paper is to study is the evolution of fold patterns that emerge out of randomly perturbed layers for different loading conditions. In order to be able to statistically analyze such systems where many folds interact, a large number of folds is required and consequently large numerical resolutions (in the order of 100’000’000 unknowns). We describe the developed numerical model and analyze the fold patterns using differential geometry. The obtained results indicate that the (Gaussian curvature based) aspect ratio of folds in map view may be used to infer the relative strength of the two principal in-plane loads.