Mikito Furuichi

Development of a Stokes flow solver robust to large viscosity jumps using a Schur complement approach with mixed precision arithmetic

We develop an iterative solution technique for solving Stokes flow problems with smooth and discontinuous viscosity structures using a three dimensional, staggered grid finite difference discretization. Two preconditioned iterative methodologies are applied to the saddle point arising from the discrete Stokes problem. They consist of a velocity-pressure coupled approach (FC) and a decoupled, Schur complement approach (SC). Within both of these methods, we utilize either the scaled BFBt, or an identity matrix scaled by the local cell viscosity (LV) to define a preconditioner for the Schur complement. Additionally, we propose to use a mixed precision Krylov kernel to improve the convergence by reducing round-off error. In this approach, standard double precision is used during the application of the preconditioner, whilst higher precision arithmetic is used to define the matrix vector product, dot products and norms required by the Krylov method. In our Krylov kernel, we utilize quad precision arithmetic which is emulated via the double-double precision method. We consider several simplified geodynamic problems with a viscosity contrast to demonstrate the robustness and scalability of our solution methods. Through a careful choice of stopping conditions, we are able to quantitatively compare the residuals between the SC and FC approaches. We examine the trade-off relationship between the number of outer iterations required for convergence, and the computational cost per iteration, for the each solution methods. We find that it is advantageous to use the FC approach utilizing relaxed tolerances for solution of the sub-problems, combined with the LV preconditioner. We also observed that in general, the SC approach is more robust than FC and that BFBt is more robust than LV when used in our numerical experimental. In addition, our mixed precision method produces improved convergence rates of Arnoldi type Krylov subspace methods without a drastic increasing the computational time. The usage of a high precision Krylov kernel is found to be useful for the solver associated with the velocity sub-problem.

Computational performance of a smoothed particle hydrodynamics simulation for shared-memory parallel computing

Abstract The computational performance of a smoothed particle hydrodynamics (SPH) simulation is investigated for three types of current shared-memory parallel computer devices: many integrated core (MIC) processors, graphics processing units (GPUs), and multi-core CPUs. We are especially interested in efficient shared-memory allocation methods for each chipset, because the efficient data access patterns differ between compute unified device architecture (CUDA) programming for GPUs and OpenMP programming for MIC processors and multi-core CPUs. We first introduce several parallel implementation techniques for the SPH code, and then examine these on our target computer architectures to determine the most effective algorithms for each processor unit. In addition, we evaluate the effective computing performance and power efficiency of the SPH simulation on each architecture, as these are critical metrics for overall performance in a multi-device environment. In our benchmark test, the GPU is found to produce the best arithmetic performance as a standalone device unit, and gives the most efficient power consumption. The multi-core CPU obtains the most effective computing performance. The computational speed of the MIC processor on Xeon Phi approached that of two Xeon CPUs. This indicates that using MICs is an attractive choice for existing SPH codes on multi-core CPUs parallelized by OpenMP, as it gains computational acceleration without the need for significant changes to the source code.

Three-dimensional Eulerian method for large deformation of viscoelastic fluid: Toward plate-mantle simulation

Toward the unified simulation of the large deformation of a rigid viscoelastic material (plate) and the convection of a viscous fluid (mantle), an Eulerian scheme with a semi-Lagrangian method is developed. The scheme adopts the CIP-CSLR method for advection terms of staggered grid system in three dimensions. The positive transported profile of a positive quantity is assured by flux corrections in the dimensional splitting method. The Jaumann co-rotational effect of the stress tensor is also integrated into the semi-Lagrangian treatment. This co-rotated semi-Lagrangian method is combined with an exponential time differencing method in the time development of the Maxwell constitutive model. The large time step comparable to, or larger than, the Maxwell relaxation time is successfully realized. Validation tests are performed for the three-dimensional Rayleigh-Taylor instability of a viscoelastic material with jump discontinuity of the mass density and other material properties.

MovieMaker: a parallel movie-making software for large-scale simulations

We have developed a parallel rendering software for scientific visualization of large-scale, three-dimensional, time development simulations. The goal of this software, MovieMaker, is to generate a movie, or a series of visualization images from one terabyte-scale data set within one night (or less than 12 h). The isorontouring, volume rendering, and streamlines are implemented. MovieMaker is a parallel program for the shared memory architecture with dynamic load balancing and overlapped disk I/O.

Implicit solution of the material transport in Stokes flow simulation: Toward thermal convection simulation surrounded by free surface

Abstract We present implicit time integration schemes suitable for modeling free surface Stokes flow dynamics with marker in cell (MIC) based spatial discretization. Our target is for example thermal convection surrounded by deformable surface boundaries to simulate the long term planetary formation process. The numerical system becomes stiff when the dynamical balancing time scale for the increasing/decreasing load by surface deformation is very short compared with the time scale associated with thermal convection. Any explicit time integration scheme will require very small time steps; otherwise, serious numerical oscillation (spurious solutions) will occur. The implicit time integration scheme possesses a wider stability region than the explicit method; therefore, it is suitable for stiff problems. To investigate an efficient solution method for the stiff Stokes flow system, we apply first (backward Euler (BE)) and second order (trapezoidal method (TR) and trapezoidal rule—backward difference formula (TR-BDF2)) accurate implicit methods for the MIC solution scheme. The introduction of implicit time integration schemes results in nonlinear systems of equations. We utilize a Jacobian free Newton Krylov (JFNK) based Newton framework to solve the resulting nonlinear equations. In this work we also investigate two efficient implicit solution strategies to reduce the computational cost when solving stiff nonlinear systems. The two methods differ in how the advective term in the material transport evolution equation is treated. We refer to the method that employs Lagrangian update as “fully implicit” (Imp), whilst the method that employs Eulerian update is referred to as “semi-implicit” (SImp). Using a finite difference (FD) method, we have performed a series of numerical experiments which clarify the accuracy of solutions and trade-off between the computational cost associated with the nonlinear solver and time step size. In comparison with the general explicit Euler method, the second order accurate Imp methods reduce total computational cost successfully through the utilization of a large time step without sacrificing accuracy and stability. Moreover, the proposed SImp method is effective in reducing the computational cost associated with evaluating the nonlinear residual while obtaining a solution similar to the Imp method.

Robust coupled fluid‐particle simulation scheme in Stokes‐flow regime: Toward the geodynamic simulation including granular media

We present a simulation scheme for solving high-viscosity fluid and particle dynamics in a coupled computational fluid dynamics and discrete element method (CFD-DEM) framework. This simulation scheme is intended to be used for geodynamical magmatic studies such as crystal settling at the melting roof of a magma chamber. The high-viscosity fluid is treated by the Stokes-flow approximation, where the fluid interacts with particles via the drag force in a cell-averaged manner. The particles are tracked with contact forces by DEM. To efficiently solve such Stokes-DEM coupled equations, we propose two key techniques. One is formulation of particle motion without the inertial term, allowing a larger time step at higher viscosities. The other is a semi-implicit treatment of the cell-averaged particle velocity in the fluid equation to stabilize the calculation. We simulate the settling particles in strongly viscous fluids in three dimensions and compare the results with the experimental and theoretical results. Our solution strategy is found to be robust and successfully captures the collective behavior of the particles. The simulation method presented here will be useful in various fields interested in long-term dynamics of high-viscosity granular media.

Iterative load-balancing method with multigrid level relaxation for particle simulation with short-range interactions

Abstract We developed dynamic load-balancing algorithms for Particle Simulation Methods (PSM) involving short-range interactions, such as Smoothed Particle Hydrodynamics (SPH), Moving Particle Semi-implicit method (MPS), and Discrete Element method (DEM). These are needed to handle billions of particles modeled in large distributed-memory computer systems. Our method utilizes flexible orthogonal domain decomposition, allowing the sub-domain boundaries in the column to be different for each row. The imbalances in the execution time between parallel logical processes are treated as a nonlinear residual. Load-balancing is achieved by minimizing the residual within the framework of an iterative nonlinear solver, combined with a multigrid technique in the local smoother. Our iterative method is suitable for adjusting the sub-domain frequently by monitoring the performance of each computational process because it is computationally cheaper in terms of communication and memory costs than non-iterative methods. Numerical tests demonstrated the ability of our approach to handle workload imbalances arising from a non-uniform particle distribution, differences in particle types, or heterogeneous computer architecture which was difficult with previously proposed methods. We analyzed the parallel efficiency and scalability of our method using Earth simulator and K-computer supercomputer systems.

Arcuate stress state in accretionary prisms from real-scale numerical sandbox experiments

The stress states in accretionary prisms are important for understanding the building and releasing of seismic energy. Numerous researchers have conducted sandbox experiments as a scaled physical analog model to understand the formation of accretionary prisms. However, measuring stress states in laboratory sandbox experiments is still practically infeasible. Here we performed real-scale numerical sandbox experiments using the discrete element method to understand the 3D stress state in the accretionary prism. Despite the nearly uniform initial conditions, macro-scale undulations of faults, which are similar to those observed in the trenches of an accretionary prism, appear. We reveal that these undulations are caused by the formation of stress arches. We show that the mechanism behind the arch formation is the discontinuous change in the stress orientation during the rearrangement of the stress chain. Furthermore, analyses demonstrate that the in-situ stress orientation from borehole data can be a signal of either the regional direction of plate convergence or the local stress orientation associated with the stress arch. The results may greatly enhance the outcome of long term monitoring in areas, such as the Nankai Trough.

Visualization analysis of large-scale three-dimensional scalar data of ocean simulation

We visualize the data set of the hind-cast experiment of OFES (Sasaki et al., 2004). The horizontal resolution and the number of vertical levels are 0.1 degree and 54, respectively. Total data size of simulated sea temperature is 3600 × 1500 × 54 Grids × 4 Bytes × 6575 Snapshots = 7.7 Tera Bytes. In our visualization results, all of the simulated events are displayed by ray-casting volume rendering in the global view like Fig. 1, which represents the IOD and El Niño simultaneously without losing the resolution of original out put data. These sets of the pictures through 54 years must be helpful for understanding the combined influences of the events.