A fully implicit model for simulating dynamo action in a Cartesian domain

We present a fully implicit numerical method to solve the incompressible MHD equations in a strongly rotating Cartesian domain. The equations are solved in a primitive variable formulation using a finite volume discretization. In order to use massively parallel computers, we applied a domain decomposition approach in space. The performance of this model is compared with an earlier model, which treated the convective terms of the equations in an explicit manner. Our results indicate that although the fully implicit method needs about three times the memory of the implicit–explicit method, it is superior in terms of computational efficiency. As an application of this model, we investigated the influence of the Prandtl number in the range of 0.01–1000 on the dynamics of the dynamo.     

Efficient fully-coupled solution techniques for two-phase flow in porous media: Parallel multigrid solution and large scale computations

This paper is concerned with the fast resolution of nonlinear and linear algebraic equations arising from a fully implicit finite volume discretization of two-phase flow in porous media. We employ a Newton-multigrid algorithm on unstructured meshes in two and three space dimensions. The discretized operator is used for the coarse grid systems in the multigrid method. Problems with discontinuous coefficients are avoided by using a newly truncated restriction operator and an outer Krylov-space method. We show an optimal order of convergence for a wide range of two-phase flow problems including heterogeneous media and vanishing capillary pressure in an experimental way. Furthermore, we present a data parallel implementation of the algorithm with speedup results.     

An implicit wave equation model for the shallow water equations

An implicit solution procedure for the wave equation form of the shallow water equations is presented. Efficiency is achieved through a Taylor expansion procedure applied to a time-varying matrix. This procedure allows matrix decompositions to be replaced by back substitutions. Isoparametric quadratic Lagrangian finite elements are employed for the spatial discretization. The Taylor expansion method is compared to different implicit and explicit solution procedures in an application to the southern part of the North Sea.     

Effects of chemical reactions on density-dependent fluid flow: on the numerical formulation and the development of instabilities

A three-dimensional, reactive numerical flow model is developed that couples chemical reactions with density-dependent mass transport and fluid flow. The model includes equilibrium reactions for the aqueous species, kinetic reactions between the solid and aqueous phases, and full coupling of porosity and permeability changes that result from precipitation and dissolution reactions in porous media. A one-step, global implicit approach is used to solve the coupled flow, transport and reaction equations with a fully implicit upstream-weighted control volume discretization. The Newton–Raphson method is applied to the discretized non-linear equations and a block ILU-preconditioned CGSTAB method is used to solve the resulting Jacobian matrix equations. This approach permits the solution of the complete set of governing equations for both concentration and pressure simultaneously affected by chemical and physical processes. A series of chemical transport simulations are conducted to investigate coupled processes of reactive chemical transport and density-dependent flow and their subsequent impact on the development of preferential flow paths in porous media. The coupled effects of the processes driving flow and the chemical reactions occurring during solute transport is studied using a carbonate system in fully saturated porous media. Results demonstrate that instability development is sensitive to the initial perturbation caused by density differences between the solute plume and the ambient groundwater. If the initial perturbation is large, then it acts as a “trigger” in the flow system that causes instabilities to develop in a planar reaction front. When permeability changes occur due to dissolution reactions occurring in the porous media, a reactive feedback loop is created by calcite dissolution and the mixed convective transport of the system. Although the feedback loop does not have a significant impact on plume shape, complex concentration distributions develop as a result of the instabilities generated in the flow system.     

Hybrid upwind discretization of nonlinear two-phase flow with gravity

Multiphase flow in porous media is described by coupled nonlinear mass conservation laws. For immiscible Darcy flow of multiple fluid phases, whereby capillary effects are negligible, the transport equations in the presence of viscous and buoyancy forces are highly nonlinear and hyperbolic. Numerical simulation of multiphase flow processes in heterogeneous formations requires the development of discretization and solution schemes that are able to handle the complex nonlinear dynamics, especially of the saturation evolution, in a reliable and computationally efficient manner. In reservoir simulation practice, single-point upwinding of the flux across an interface between two control volumes (cells) is performed for each fluid phase, whereby the upstream direction is based on the gradient of the phase-potential (pressure plus gravity head). This upwinding scheme, which we refer to as Phase-Potential Upwinding (PPU), is combined with implicit (backward-Euler) time discretization to obtain a Fully Implicit Method (FIM). Even though FIM suffers from numerical dispersion effects, it is widely used in practice. This is because of its unconditional stability and because it yields conservative, monotone numerical solutions. However, FIM is not unconditionally convergent. The convergence difficulties are particularly pronounced when the different immiscible fluid phases switch between co-current and counter-current states as a function of time, or (Newton) iteration. Whether the multiphase flow across an interface (between two control-volumes) is co-current, or counter-current, depends on the local balance between the viscous and buoyancy forces, and how the balance evolves in time. The sensitivity of PPU to small changes in the (local) pressure distribution exacerbates the problem. The common strategy to deal with these difficulties is to cut the timestep and try again. Here, we propose a Hybrid-Upwinding (HU) scheme for the phase fluxes, then HU is combined with implicit time discretization to yield a fully implicit method. In the HU scheme, the phase flux is divided into two parts based on the driving force. The viscous-driven and buoyancy-driven phase fluxes are upwinded differently. Specifically, the viscous flux, which is always co-current, is upwinded based on the direction of the total-velocity. The buoyancy-driven flux across an interface is always counter-current and is upwinded such that the heavier fluid goes downward and the lighter fluid goes upward. We analyze the properties of the Implicit Hybrid Upwinding (IHU) scheme. It is shown that IHU is locally conservative and produces monotone, physically-consistent numerical solutions. The IHU solutions show numerical diffusion levels that are slightly higher than those for standard FIM (i.e., implicit PPU). The primary advantage of the IHU scheme is that the numerical overall-flux of a fluid phase remains continuous and differentiable as the flow regime changes between co-current and counter-current conditions. This is in contrast to the standard phase-potential upwinding scheme, in which the overall fractional-flow (flux) function is non-differentiable across the boundary between co-current and counter-current flows.     

Flux-Corrected Transport with MT3DMS for Positive Solution of Transport with Full-Tensor Dispersion

Solute transport is usually modeled by the advection-dispersion-reaction equation. In the standard approach, mechanical dispersion is a tensor with principal directions parallel and perpendicular to the flow vector. Since realistic scenarios include nonuniform and unsteady flow fields, the governing equation has full tensor mechanical dispersion. When conventional grid-based numerical methods are used, approximation of the cross terms arising from the off-diagonal terms cause nonphysical solution with oscillations. As an example, for the common scenario of contaminant input into a domain with zero initial concentration, the cross-dispersion terms can result in negative concentrations that can wreak havoc in reactive transport applications. To address this issue, we use the well-known flux-corrected-transport (FCT) technique for a standard finite volume method. Although FCT has most often been used to eliminate oscillations resulting from discretization of the advection term for explicit time stepping, we show that it can be adapted for full-tensor dispersion and implicit time stepping. Unlike other approaches based on new discretization techniques (e.g., mimetic finite difference, nonlinear finite volume), FCT has the advantage of being flexible and widely applicable. Implementation of FCT requires solving an additional system of equations at each time step, using a modified “low order” matrix and a modified right-hand-side vector. To demonstrate the flexibility of FCT, we have modified the well-known and widely used groundwater solute transport simulator, MT3DMS. We apply the new simulator, MT3DMS-FCT, to several benchmark problems that suffer from negative concentrations when using MT3DMS. The new results are mass conservative and strictly nonnegative.     

Applying Higher Order DIRK Time Steps to the "Modified Picard" Method

The “modified Picard” iteration method, which offers global mass conservation, can also be described as a form of Newton's iteration with lagged nonlinear coefficients. It converges to a time step with first-order discretization error. This paper applies second- and third-order diagonally implicit Runge Kutta (DIRK) time steps to the modified Picard method in one example. It demonstrates improvements over the first-order time step in rms error and error-times-effort model quality by factors ranging from two to over two orders of magnitude, showing that the “modified Picard” and DIRK methods are compatible.     

Perfectly matched layer boundary conditions for frequency-domain acoustic wave simulation in the mesh-free discretization

Mesh-free discretization, flexibly distributing nodes without computationally expensive meshing process, is able to deal with staircase problem, oversampling and undersampling problems and saves plenty of nodes through distributing nodes suitably with respect to irregular boundaries and model parameters. However, the time-domain mesh-free discretization usually exhibits poorer stability than that in regular grid discretization. In order to reach unconditional stability and easy implementation in parallel computing, we develop the frequency-domain finite-difference method in a mesh-free discretization, incorporated with two perfectly matched layer boundary conditions. Furthermore, to maintain the flexibility of mesh-free discretization, the nodes are still irregularly distributed in the absorbing zone, which complicates the situation of artificial boundary reflections. In this paper, we implement frequency-domain acoustic wave modelling in a mesh-free system. First, we present the perfectly matched layer boundary condition to suppress spurious reflections. Moreover, we develop the complex frequency shifted–perfectly matched layer boundary condition to improve the attenuation of grazing waves. In addition, we employ the radial-basis-function-generated finite difference method in the mesh-free discretization to calculate spatial derivatives. The numerical experiment on a rectangle homogeneous model shows the effectiveness of the perfectly matched layer boundary condition and the complex frequency shifted–perfectly matched layer boundary condition, and the latter one is better than the former one when absorbing large angle incident waves. The experiment on the Marmousi model suggests that the complex frequency shifted–perfectly matched layer boundary condition works well for complicated models.     

三维时间域航空电磁有理Krylov正演研究

常规的三维时间域航空电磁模拟通常采用隐式步长方法进行时间离散,需要几次矩阵分解和上百次右端源项回带,计算效率较低.为了提高正演计算效率,本文提出使用有理Krylov方法求解时间域电场扩散方程.首先使用非结构四面体网格进行空间离散,采用Nédélec矢量基函数近似四面体单元内的电场;然后基于有限元离散给出矩阵指数和矢量乘积表示的电场显式解;最后采用有理Arnoldi算法构造Krylov子空间内的正交基函数并进一步求解矩阵指数与矢量的乘积,直接得到任意时刻的电场解向量,避免步长离散过程.此外,本文还提出一种指数加权偏移参数优化方法,使得有理Arnoldi近似在瞬变衰减晚期具备更高的精度,从而降低Krylov子空间阶数并提高计算效率.通过和层状模型解析解的对比验证了有理Krylov方法的精度.针对三维异常体模型使用全局网格和局部网格剖分并和其他数值方法比较,进一步说明了有理Krylov方法的有效性.     

Iterative implicit integration procedure for hybrid simulation of large nonlinear structures

A fully implicit iterative integration procedure is presented for local and geographically distributed hybrid simulation of the seismic response of complex structural systems with distributed nonlinear behavior. The purpose of this procedure is to seamlessly incorporate experimental elements in simulations using existing fully implicit integration algorithms designed for pure numerical simulations. The difficulties of implementing implicit integrators in a hybrid simulation are addressed at the element level by introducing a safe iteration strategy and using an efficient procedure for online estimation of the experimental tangent stiffness matrix. In order to avoid physical application of iterative displacements, the required experimental restoring force at each iteration is estimated from polynomial curve fitting of recent experimental measurements. The experimental tangent stiffness matrix is estimated by using readily available experimental measurements and by a classical diagonalization approach that reduces the number of unknowns in the matrix. Numerical and hybrid simulations are used to demonstrate that the proposed procedure provides an efficient method for implementation of fully implicit numerical integration in hybrid simulations of complex nonlinear structures. The hybrid simulations presented include distributed nonlinear behavior in both the numerical and experimental substructures. Copyright © 2010 John Wiley & Sons, Ltd.     

钢筋混凝土框架结构计算机仿真与并行计算

介绍了钢筋混凝土框架结构计算机仿真和并行计算的研究现状。数值仿真主要采用有限元、离散元等数学物理模型,而可视化技术及图形仿真基于图形学和图像处理技术,是驾驭计算过程及理解大体积数据的唯一有效途径。有限元并行计算有SBS、EBE两种策略,非线性求解有直接与迭代解法,动力时程分析有显式、隐式和精细时程积分法。最后介绍了并行计算在钢筋混凝土结构分析中的应用。     

A parallel,fully coupled,fully implicit solution to reactive transport in porous media using the preconditioned Jacobian-Free Newton-Krylov Method

Modeling large multicomponent reactive transport systems in porous media is particularly challenging when the governing partial differential algebraic equations (PDAEs) are highly nonlinear and tightly coupled due to complex nonlinear reactions and strong solution-media interactions. Here we present a preconditioned Jacobian-Free Newton-Krylov (JFNK) solution approach to solve the governing PDAEs in a fully coupled and fully implicit manner. A well-known advantage of the JFNK method is that it does not require explicitly computing and storing the Jacobian matrix during Newton nonlinear iterations. Our approach further enhances the JFNK method by utilizing physics-based, block preconditioning and a multigrid algorithm for efficient inversion of the preconditioner. This preconditioning strategy accounts for self- and optionally, cross-coupling between primary variables using diagonal and off-diagonal blocks of an approximate Jacobian, respectively. Numerical results are presented demonstrating the efficiency and massive scalability of the solution strategy for reactive transport problems involving strong solution-mineral interactions and fast kinetics. We found that the physics-based, block preconditioner significantly decreases the number of linear iterations, directly reducing computational cost; and the strongly scalable algebraic multigrid algorithm for approximate inversion of the preconditioner leads to excellent parallel scaling performance.     

A Douglas-Jones Predictor-Corrector Program for Simulating One-Dimensional Unsaturated Flow in Soil

Abstract. A fully documented program to represent one-dimensional unsaturated flow in soil is described. The program is based on a Douglas-Jones finite-difference implicit method to solve the Richards equation. An implicit linearization scheme is used to estimate the hydraulic conductivity and specific moisture capacity functions. Predicted moisture content profiles compared with two Galerkin finite-element solutions and field observations on a Panoche clay loam soil show very good agreement.     

Effects of chemical reactions on iterative methods for implicit time stepping

We compare the performance of the fully coupled Newton–Raphson method with the sequential iteration approach (SIA) for solving the implicit time stepping equations of reactive transport modeling. We formulate the implicit time stepping equations for a demonstration model that incorporates homogeneous equilibrium reactions, i.e. carbonate hydrolysis, and a heterogeneous equilibrium reaction, i.e. the dissolution/precipitation of calcite. The demonstration model uses a coupled pair of mixing cells as a simplified form of transport. The effects of the homogeneous and heterogeneous reactions on the iterative methods are demonstrated.     

A terrain-following grid transform and preconditioner for parallel,large-scale,integrated hydrologic modeling

A terrain-following grid formulation (TFG) is presented for simulation of coupled variably-saturated subsurface and surface water flow. The TFG is introduced into the integrated hydrologic model, ParFlow, which uses an implicit, Newton Krylov solution technique. The analytical Jacobian is also formulated and presented and both the diagonal and non-symmetric terms are used to precondition the Krylov linear system. The new formulation is verified against an orthogonal stencil and is shown to provide increased accuracy at lower lateral spatial discretization for hillslope simulations. Using TFG, efficient scaling to a large number of processors (16,384) and a large domain size (8.1 Billion unknowns) is shown. This demonstrates the applicability of this formulation to high-resolution, large-spatial extent hydrology applications where topographic effects are important. Furthermore, cases where the analytical Jacobian is used for the Newton iteration and as a non-symmetric preconditioner for the linear system are shown to have faster computation times and better scaling. This demonstrates the importance of solver efficiency in parallel scaling through the use of an appropriate preconditioner.     

基于MPICH环境下的复杂重力异常体并行正演模拟

任意多面体重力异常正演公式常用于解决复杂几何形态地质体的正演问题。 本文以均匀物性多面体重力异常正演公式为基础, 应用有限元技术中的网格离散化思想, 以任意四面体为基本单元, 通过并行计算技术在MPICH环境下实现了任意连续空间物性分布复杂异常体网格模型的重力异常正演模拟, 通过并行处理可以有效加速正演计算速度。 本文研究结果对于联合重力异常场正演建模和开展复杂模型网格的重力场计算有一定参考意义。     

A parallel laboratory for simulation and visualization of seismic wavefields

The use of parallel computers makes simulation of elastic waves feasible throughout large structures by means of recent advances in domain decomposition methods. We introduce a competitive parallel algorithm for the propagation of elastic waves in complex heterogeneous media using finite-element discretization. This parallel method, called the multiblock method , performs more efficiently than classical domain decomposition techniques based on substructuration, such as the Schur complement technique. It reduces considerably the amount of communication amongst processors because the interface problem between subdomains is solved by taking advantage of Huygens' principle for wave propagation. We provide some numerical examples and detailed studies on the efficiency and performance of the algorithm, proving that it is competitive and less costly, from the computational viewpoint, than algorithms based on the Schur technique.     

A comparison of a semi-implicit with an explicit scheme in a three-dimensional hydrodynamic model

As the capacity of computers increases, the size and resolution of numerical models can be increased. In tidal models, however, using an explicit scheme together with high spatial resolution results in an unreasonably small time-step demanded by the CFL condition for stability. This condition is usually the most restrictive and applies to the propagation of the gravity wave.A fully nonlinear three-dimensional model has been developed, using the Galerkin method in the vertical dimension, in which the gravity wave terms are treated by an alternating-direction implicit scheme, while the friction, viscous and advective terms are treated explicitly. This permits a stable solution with a longer time-step than that required in an equivalent explicit scheme, while not needing as much computational effort as a fully implicit scheme. This semi-implicit model is compared with an explicit model in terms of efficiency, accuracy, and stability. Tidal and wind-driven flows and free oscillations in a rectangular sea model are examined, using various boundary conditions. The semi-implicit scheme takes approximately 1.1 times as long to run (per time-step) on the CRAY-1 computer as the explicit scheme.     

Modeling Surface Water-Groundwater Interaction with MODFLOW: Some Considerations

The accuracy with which MODFLOW simulates surface water-groundwater interaction is examined for connected and disconnected losing streams. We compare the effect of different vertical and horizontal discretization within MODFLOW and also compare MODFLOW simulations with those produced by HydroGeoSphere. HydroGeoSphere is able to simulate both saturated and unsaturated flow, as well as surface water, groundwater and the full coupling between them in a physical way, and so is used as a reference code to quantify the influence of some of the simplifying assumptions of MODFLOW. In particular, we show that (1) the inability to simulate negative pressures beneath disconnected streams in MODFLOW results in an underestimation of the infiltration flux; (2) a river in MODFLOW is either fully connected or fully disconnected, while in reality transitional stages between the two flow regimes exist; (3) limitations in the horizontal discretization of the river can cause a mismatch between river width and cell width, resulting in an error in the water table position under the river; and (4) because coarse vertical discretization of the aquifer is often used to avoid the drying out of cells, this may result in an error in simulating the height of the groundwater mound. Conditions under which these errors are significant are investigated.     

Proper orthogonal decomposition reduced model for mass transport in heterogenous media

Numerical models with fine discretization normally demand large computational time and space, which lead to computational burden for state estimations or model parameter inversion calculation. This article presented a reduced implicit finite difference scheme that based on proper orthogonal decomposition (POD) for two-dimensional transient mass transport in heterogeneous media. The reduction of the original full model was achieved by projecting the high-dimension full model to a low-dimension space created by POD bases, and the bases are derived from the snapshots generated from the model solutions of the forward simulations. The POD bases were extracted from the ensemble of snapshots by singular value decomposition. The dimension of the Jacobian matrix was then reduced after Galerkin projection. Thus, the reduced model can accurately reproduce and predict the original model’s transport process with significantly decreased computational time. This scheme is practicable with easy implementation of the partial differential equations. The POD method is illustrated and validated through synthetic cases with various heterogeneous permeability field scenarios. The accuracy and efficiency of the reduced model are determined by the optimal selection of the snapshots and POD bases.