A parallel multiphase flow code for the 3D simulation of explosive volcanic eruptions

A new parallel code for the simulation of the transient, 3D dispersal of volcanic particles in the atmosphere is presented. The model equations, describing the multiphase flow dynamics of gas and solid pyroclasts ejected from the volcanic vent during explosive eruptions, are solved by a finite-volume discretization scheme and a pressure-based iterative non-linear solver suited to compressible multiphase flows. The solution of the multiphase equation set is computationally so demanding that the simulation of the transient 3D dynamics of eruptive columns would not be cost-effective on a single workstation. The new code has been parallelized by adopting an ad hoc domain partitioning scheme that enforces the load balancing in the presence of a large number of topographic blocking-cells. An optimized communication layer has been built over the Message-Passing Interface. It is shown that the present code has a remarkable efficiency on several high-performance platforms and makes it possible, for the first time, to simulate fully 3D eruptive scenarios on realistic volcano topography.     

Structure topology optimization: fully coupled level set method via FEMLAB

This paper presents a procedure which can easily implement the 2D compliance minimization structure topology optimization by the level set method using the FEMLAB package. Instead of a finite difference solver for the level set equation, as is usually the case, a finite element solver for the reaction–diffusion equation is used to evolve the material boundaries. All of the optimization procedures are implemented in a user-friendly manner. A FEMLAB code can be downloaded from the homepage www.imtek.de/simulation and is free for educational purposes.     

A block LU-SGS implicit unsteady incompressible flow solver on hybrid dynamic grids for 2D external bio-fluid simulations

A hybrid dynamic grid generation technique for two-dimensional (2D) morphing bodies and a block lower-upper symmetric Gauss-Seidel (BLU-SGS) implicit dual-time-stepping method for unsteady incompressible flows are presented for external bio-fluid simulations. To discretize the complicated computational domain around 2D morphing configurations such as fishes and insect/bird wings, the initial grids are generated by a hybrid grid strategy firstly. Body-fitted quadrilateral (quad) grids are generated first near solid bodies. An adaptive Cartesian mesh is then generated to cover the entire computational domain. Cartesian cells which overlap the quad grids are removed from the computational domain, and a gap is produced between the quad grids and the adaptive Cartesian grid. Finally triangular grids are used to fill this gap. During the unsteady movement of morphing bodies, the dynamic grids are generated by a coupling strategy of the interpolation method based on ‘Delaunay graph’ and local remeshing technique. With the motion of moving/morphing bodies, the grids are deformed according to the motion of morphing body boundaries firstly with the interpolation strategy based on ‘Delaunay graph’ proposed by Liu and Qin. Then the quality of deformed grids is checked. If the grids become too skewed, or even intersect each other, the grids are regenerated locally. After the local remeshing, the flow solution is interpolated from the old to the new grid. Based on the hybrid dynamic grid technique, an efficient implicit finite volume solver is set up also to solve the unsteady incompressible flows for external bio-fluid dynamics. The fully implicit equation is solved using a dual-time-stepping approach, coupling with the artificial compressibility method (ACM) for incompressible flows. In order to accelerate the convergence history in each sub-iteration, a block lower-upper symmetric Gauss-Seidel implicit method is introduced also into the solver. The hybrid dynamic grid generator is tested by a group of cases of morphing bodies, while the implicit unsteady solver is validated by typical unsteady incompressible flow case, and the results demonstrate the accuracy and efficiency of present solver. Finally, some applications for fish swimming and insect wing flapping are carried out to demonstrate the ability for 2D external bio-fluid simulations.     

A scalable parallel Poisson solver for three-dimensional problems with one periodic direction

A code for the direct numerical simulation (DNS) of incompressible flows with one periodic direction has been developed. It provides a fairly good performance on both Beowulf clusters and supercomputers. Since the code is fully explicit, from a parallel point-of-view, the main bottleneck is the Poisson equation. To solve it, a Fourier diagonalization is applied in the periodic direction to decompose the original 3D system into a set of mutually independent 2D systems. Then, different strategies can be used to solved them. In the previous version of the code, that was conceived for low-cost PC clusters with poor network performance, a Direct Schur-complement Decomposition (DSD) algorithm was used to solve them. Such a method, that is very efficient for PC clusters, cannot be used with an arbitrarily large number of processors and mesh sizes, mainly due to the RAM memory requirements. To do so, a new version of the solver is presented in this paper. It is based on the DSD algorithm that is used as a preconditioner for a Conjugate Gradient method. Numerical experiments showing the scalability and the flexibility of the method on both the MareNostrum supercomputer and a PC cluster with a conventional 100 Mbits/s network are presented and discussed. Finally, illustrative DNS results of an air-filled differentially heated cavity at Ra = 10¹¹ are also presented.     

Parallelization and optimization of electrostatic Particle-in-Cell/Monte-Carlo Coupled codes as applied to RF discharges

We have developed a full paralleled 2D electrostatic Particle-in-Cell/Monte-Carlo Coupled (PIC-MCC) code for capacitively coupled plasma (CCP) simulations. In this code, we distributed the grid between processors along radial direction, and Poisson equation is solved accordingly paralleled. We applied a couple of numerical accelerating technologies: paralleled fast Poisson solver, assembler pushing code, particle sorting and so on. Theoretical analysis and numerical benchmark showed that this parallel framework had good efficiency and scalability. The framework of the code and the optimization technologies and algorithms are discussed, benchmarks and simulation results are also shown.     

Density based topology optimization of turbulent flow heat transfer systems

The focus of this article is on topology optimization of heat sinks with turbulent forced convection. The goal is to demonstrate the extendibility, and the scalability of a previously developed fluid solver to coupled multi-physics and large 3D problems. The gradients of the objective and the constraints are obtained with the help of automatic differentiation applied on the discrete system without any simplifying assumptions. Thus, as demonstrated in earlier works of the authors, the sensitivities are exact to machine precision. The framework is applied to the optimization of 2D and 3D problems. Comparison between the simplified 2D setup and the full 3D optimized results is provided. A comparative study is also provided between designs optimized for laminar and turbulent flows. The comparisons highlight the importance and the benefits of full 3D optimization and including turbulence modeling in the optimization process, while also demonstrating extension of the methodology to include coupling of heat transfer with turbulent flows.     

An open source, parallel DSMC code for rarefied gas flows in arbitrary geometries

This paper presents the results of validation of an open source Direct Simulation Monte Carlo (DSMC) code for general application to rarefied gas flows. The new DSMC code, called dsmcFoam, has been written within the framework of the open source C++ CFD toolbox OpenFOAM. The main features of dsmcFoam code include the capability to perform both steady and transient solutions, to model arbitrary 2D/3D geometries, and unlimited parallel processing. Test cases have been selected to cover a wide range of benchmark examples from 1D to 3D. These include relaxation to equilibrium, 2D flow over a flat plate and a cylinder, and 3D supersonic flows over complex geometries. In all cases, dsmcFoam shows very good agreement with data provided by both analytical solutions and other contemporary DSMC codes.     

Parallel multigrid finite volume computation of three-dimensional thermal convection

A parallel implementation of the finite volume method for three-dimensional, time-dependent, thermal convective flows is presented. The algebraic equations resulting from the finite volume discretization, including a pressure equation which consumes most of the computation time, are solved by a parallel multigrid method. A flexible parallel code has been implemented on the Intel Paragon, the Cray T3D, and the IBM SP2 by using domain decomposition techniques and the MPI communication software. The code can use 1D, 2D, or 3D partitions as required by different geometries, and is easily ported to other parallel systems. Numerical solutions for air (Prandtl number Pr = 0.733) with various Rayleigh numbers up to 10⁷ are discussed.     

Extension of a hybrid thermal LBE scheme for large-eddy simulations of turbulent convective flows

Following the work of Lallemand and Luo [Lallemand P, Luo L-S. Theory of the lattice Boltzmann method: acoustic and thermal properties in two and three dimensions. Phys Rev E 2003;68:036706] we validate, apply and extend the hybrid thermal lattice Boltzmann scheme (HTLBE) by a large-eddy approach to simulate turbulent convective flows. For the mass and momentum equations, a multiple-relaxation-time LBE scheme is used while the heat equation is solved numerically by a finite difference scheme. We extend the hybrid model by a Smagorinsky subgrid scale model for both the fluid flow and the heat flux. Validation studies are presented for laminar and turbulent natural convection in a cavity at various Rayleigh numbers up to 5 × 10¹⁰ for Pr = 0.71 using a serial code in 2D and a parallel code in 3D, respectively. Correlations of the Nusselt number are discussed and compared to benchmark data. As an application we simulated forced convection in a building with inner courtyard at Re = 50 000.     

Efficient Parallelization of a Three-Dimensional Navier-Stokes Solver on MIMD Multiprocessors

We present an efficient parallelization strategy for speeding up the computation of a high-accuracy 3-dimensional serial Navier-Stokes solver that treats turbulent transonic high-Reynolds flows. The code solves the full compressible Navier-Stokes equations and is applicable to realistic large size aerodynamic configurations and as such requires huge computational resources in terms of computer memory and execution time. The solver can resolve the flow properly on relatively coarse grids. Since the serial code contains a complex infrastructure typical for industrial code (which ensures its flexibility and applicability to complex configurations), then the parallelization task is not straightforward. We get scalable implementation on massively parallel machines by maintaining efficiency at a fixed value by simultaneously increasing the number of processors and the size of the problem.The 3-D Navier-Stokes solver was implemented on three MIMD message-passing multiprocessors (a 64-processors IBM SP2, a 20-processors MOSIX, and a 64-processors Origin 2000). The same code written with PVM and MPI software packages was executed on all the above distinct computational platforms. The examples in the paper demonstrate that we can achieve efficiency of about 60% for as many as 64 processors on Origin 2000 on a full-size 3-D aerodynamic problem which is solved on realistic computational grids.     

Development of a parallel implicit solver of fluid modeling equations for gas discharges

A parallel fully implicit PETSc-based fluid modeling equations solver for simulating gas discharges is developed. Fluid modeling equations include: the neutral species continuity equation, the charged species continuity equation with drift-diffusion approximation for mass fluxes, the electron energy density equation, and Poisson's equation for electrostatic potential. Except for Poisson's equation, all model equations are discretized by the fully implicit backward Euler method as a time integrator, and finite differences with the Scharfetter–Gummel scheme for mass fluxes on the spatial domain. At each time step, the resulting large sparse algebraic nonlinear system is solved by the Newton–Krylov–Schwarz algorithm. A 2D-GEC RF discharge is used as a benchmark to validate our solver by comparing the numerical results with both the published experimental data and the theoretical prediction. The parallel performance of the solver is investigated.     

An implicit compact scheme solver for two-dimensional multicomponent flows

A 2D implicit compact scheme solver has been implemented for the vorticity-velocity formulation in the case of nonreacting, multicomponent, axisymmetric, low Mach number flows. To stabilize the discrete boundary value problem, two sets of boundary closures are introduced to couple the velocity and vorticity fields. A Newton solver is used for solving steady-state and time-dependent equations. In this solver, the Jacobian matrix is formulated and stored in component form. To solve the system of linearized equations within each iteration of Newton’s method, preconditioned Bi-CGSTAB is used in combination with a matrix-vector product computed in component form. The almost dense Jacobian matrix is approximated by a partial Jacobian. For the preconditioner equation, the partial Jacobian is approximately factored using several methods. In a detailed study of several preconditioning techniques, a promising method based on ILUT preconditioning in combination with reordering and double scaling using the MC64 algorithm by Duff and Koster is selected. To validate the implicit compact scheme solver, several nonreacting model problems have been considered. At least third order accuracy in space is recovered on nonuniform grids. A comparison of the results of the implicit compact scheme solver with the results of a traditional implicit low order solver shows an order of magnitude reduction of computer memory and time using the compact scheme solver in the case of time-dependent mixing problems.     

Parallelizaton of visual magneto-hydrodynamics code based on fluctuation distribution scheme on triangular grids

A two-dimensional (2D) magneto-hydrodynamics (MHD) code which has visualization and parallel processing capability is presented in this paper. The code utilizes a fluctuation splitting (FS) scheme that runs on structured or unstructured triangular meshes. First FS scheme which included the wave model: Model-A had been developed by Roe [Roe PL. Discrete models for the numerical analysis of time-dependent multi-dimensional gas dynamics. J Comp Phys 1986;63:458-76.] for the solutions of Euler’s equations. The first 2D-MHD wave model: MHD-A, was then developed by Balci and Aslan [Balci ?. The numerical solutions of two dimensional MHD equations by fluctuation splitting scheme on triangular meshes, Ph.D. Thesis, University of Marmara, Science-Art Faculty, Physics Dept Istanbul, Turkey; 2000; Aslan N. MHD-A: A fluctuation splitting wave model for planar magnetohydrodynamics. J Comp Phys 1999;153:437-66.] to solve MHD problems including shocks and discontinuities. It was then shown in [Balci S, Aslan N. Two dimensional MHD solver by fluctuation splitting and dual time stepping. Int J Numer Meth Fluids, in press.] that this code was capable of producing reliable results in compressible and nearly incompressible limits and under the effect of gravitational fields and that it was able to identically reduce to model-A of Roe in Euler limit with no sonic problems at rarefaction fans (Balci and Aslan, in press). An important feature of this code is its ability to run time dependent or steady problems on structured or unstructured triangular meshes that can be generated automatically by the code for specified domains. In order to use the parallel processing capability of the code, the triangular meshes are decomposed into different blocks in order to share the workload among a number of processors (here personal computers) which are connected by Ethernet. Due to the compact nature of the FS scheme, only one set of data transfer is required between neighbor processors. As it will be shown, this phenomenon results in minimum amount of communication loss and makes the scheme rather robust for parallel processing. The other important feature of the new code is its visual capability. As the code is running, colorful images of scalar quantities (density, pressure, Mach number, etc.) or vector graphics of vectoral quantities (velocity, magnetic field, etc.) can be followed on the screen. The extended code, called PV-MHDA, also allows following the trajectories of the particles in time by means of a recently included particle in cell (PIC) algorithm. Because the numerical dissipation embedded in its wave model reflects real physical viscosity and resistivity, it is able to run accurately for compressible flows (including shocks) as well as nearly incompressible flows (e.g., Kelvin-Helmholtz instability). The user-friendly visual and large-scale computation capability of the code allow the user more thorough analysis of MHD problems in two-dimensional complex domains.     

A global collisionless PIC code in magnetic coordinates

A global plasma turbulence simulation code, ORB5, is presented. It solves the gyrokinetic electrostatic equations including zonal flows in axisymmetric magnetic geometry. The present version of the code assumes a Boltzmann electron response on magnetic surfaces. It uses a Particle-In-Cell (PIC), δf scheme, 3D cubic B-splines finite elements for the field solver and several numerical noise reduction techniques. A particular feature is the use of straight-field-line magnetic coordinates and a field-aligned Fourier filtering technique that dramatically improves the performance of the code in terms of both the numerical noise reduction and the maximum time step allowed. Another feature is the capability to treat arbitrary axisymmetric ideal MHD equilibrium configurations. The code is heavily parallelized, with scalability demonstrated up to 4096 processors and 10⁹ marker particles. Various numerical convergence tests are performed. The code is validated against an analytical theory of zonal flow residual, geodesic acoustic oscillations and damping, and against other codes for a selection of linear and nonlinear tests.     

A Simple Solver for Linear Equations Containing Nonlinear Operators

This paper presents a simple equation solver. The solver finds solutions for sets of linear equations extended with several nonlinear operators, including integer division and modulus, sign extension, and bit slicing. The solver uses a new technique called {\em balancing}, which can eliminate some nonlinear operators from a set of equations before applying Gaussian elimination. The solver's principal advantages are its simplicity and its ability to handle some nonlinear operators, including nonlinear functions of more than one variable. The solver is part of an application generator that provides encoding and decoding of machine instructions based on equational specifications. The solver is presented not as pseudo code but as a literate program, which guarantees that the code shown in the paper is the same code that is actually used. Using real code exposes more detail than using pseudocode, but literate-programming techniques help manage the detail. The detail should benefit readers who want to implement their own solvers based on the techniques presented here.     

Lattice Boltzmann equation for axisymmetric thermal flows

A simple lattice Boltzmann equation (LBE) model for axisymmetric thermal flow is proposed in this paper. The flow field is solved by a quasi-two-dimensional nine-speed (D2Q9) LBE, while the temperature field is solved by another four-speed (D2Q4) LBE. The model is validated by a thermal flow in a pipe and some nontrivial thermal buoyancy-driven flows in vertical cylinders, including Rayleigh-Bénard convection, natural convection, and heat transfer of swirling flows. It is found that the numerical results agree excellently with analytical solution or other numerical results.     

A p-multigrid spectral difference method for two-dimensional unsteady incompressible Navier-Stokes equations

This paper presents the development of a 2D high-order solver with spectral difference method for unsteady incompressible Navier-Stokes equations accelerated by a p-multigrid method. This solver is designed for unstructured quadrilateral elements. Time-marching methods cannot be applied directly to incompressible flows because the governing equations are not hyperbolic. An artificial compressibility method (ACM) is employed in order to treat the inviscid fluxes using the traditional characteristics-based schemes. The viscous fluxes are computed using the averaging approach (Sun et al., 2007; Kopriva, 1998) [29] and [12]. A dual time stepping scheme is implemented to deal with physical time marching. A p-multigrid method is implemented (Liang et al., 2009) [16] in conjunction with the dual time stepping method for convergence acceleration. The incompressible SD (ISD) method added with the ACM (SD-ACM) is able to accurately simulate 2D steady and unsteady viscous flows.     

A DEM-DLM/FD method for direct numerical simulation of particulate flows: Sedimentation of polygonal isometric particles in a Newtonian fluid with collisions

An efficient direct numerical simulation method to tackle the problem of particulate flows at moderate to high concentration and finite Reynolds number is presented. Our method is built on the framework established by Glowinski and his co-workers [Glowinski R, Pan TW, Hesla TI, Joseph DD. A distributed lagrange multiplier/fictitious domain method for particulate flow. Int J Multiphase Flow 1999;25:755-94] in the sense that we use their Distributed Lagrange Multiplier/Fictitious Domain (DLM/FD) formulation and their operator-splitting idea but differs in the treatment of particle collisions. Compared to our previous works [Yu Z, Wachs A, Peysson Y. Numerical simulation of particle sedimentation in shear-thinning fluids with a fictitious domain method. J Non Newtonian Fluid Mech 2006;136:126-139; Yu Z, Shao X, Wachs A. A fictitious domain method for particulate flow with heat transfer. J Comput Phys 2006;217:424-52; Yu Z, Wachs A. A fictitious domain method for dynamic simulation of particle sedimentation in Bingham fluids. J Non Newtonian Fluid Mech 2007;145:78-91], the novelty of our present contribution relies on replacing the simple artificial repulsive force based collision model usually employed in the literature by an efficient Discrete Element Method (DEM) granular solver. The use of our DEM solver enables us to consider particles of arbitrary shape (at least convex) and to account for actual contacts, in the sense that particles actually touch each other, in contrast with the repulsive force based collision model. We validate GRIFF,¹ our numerical code, against benchmark problems and compare our predictions with those available in the literature. Results, which, to the best of our knowledge, have never been reported elsewhere, on the 2D sedimentation of isometric polygonal particles with collisions are presented.     

A 3D finite volume method for the prediction of a supercritical fluid buoyant flow in a differentially heated cavity

This paper describes a three-dimensional finite volume method for the prediction of supercritical fluid buoyant flows in heated enclosures. The space and the time accuracies of the method are checked on an exact analytical solution, and the solver is validated for several benchmark tests for natural convection inside a differentially heated cavity both in the Boussinesq and in the low Mach number approximations. Then, the influence of the linearization of the van der Waals equation of state on the global convergence is discussed. Comparisons, based on the Rayleigh number, between the supercritical fluid and the perfect gas flows in a side heated cavity, show many similarities in the thermal equilibria and considerable differences in the transient behaviours.