3D domain decomposition for violent wave‐ship interactions

A 3D Domain‐Decomposition (DD) strategy has been developed to deal with violent wave‐ship interactions involving water‐on‐deck and slamming occurrence. It couples a linear potential flow seakeeping solver with a Navier–Stokes method. The latter is applied in an inner domain where slamming, water‐on‐deck, and free surface fragmentation may occur, involving important flow nonlinearities. The field solver combines an approximated projection method with a level set technique for the free surface evolution. A hybrid strategy, combining the Eulerian level set concept to Lagrangian markers, is used to enforce more accurately the body boundary condition in case of high local curvatures. Main features of the weak and the strong coupling algorithms are described with special focus on the boundary conditions for the inner solver. Two ways of estimating the nonlinear loads by the Navier–Stokes method are investigated, on the basis of an extrapolation technique and an interpolation marching cubes algorithm, respectively. The DD is applied for the case of a freely floating patrol ship in head sea regular waves and compared against water‐on‐deck experiments in terms of flow evolution, body motions, and pressure on the hull. Improvement of the solver efficiency and accuracy is suggested. Copyright © 2013 John Wiley & Sons, Ltd.     

Newton‐Krylov method for computing the cyclic steady states of evolution problems in nonlinear mechanics

This work is focused on the Newton‐Krylov technique for computing the steady cyclic states of evolution problems in nonlinear mechanics with space‐time periodicity conditions. This kind of problems can be faced, for instance, in the modeling of a rolling tire with a periodic tread pattern, where the cyclic state satisfies “rolling” periodicity condition, including shifts both in time and space. The Newton‐Krylov method is a combination of a Newton nonlinear solver with a Krylov linear solver, looking for the initial state, which provides the space‐time periodic solution. The convergence of the Krylov iterations is proved to hold in presence of an adequate preconditioner. After preconditioning, the Newton‐Krylov method can be also considered as an observer‐controller method, correcting the transient solution of the initial value problem after each period. Using information stored while computing the residual, the Krylov solver computation time becomes negligible with respect to the residual computation time. The method has been analyzed and tested on academic applications and compared with the standard evolution (fixed point) method. Finally, it has been implemented into the Michelin industrial code, applied to a full 3D rolling tire model.     

Parallel computation of a damage localization problem using parallel multifrontal solver

In this paper, nonlinear parallel structural analyses are performed using a distributed memory sparse direct multifrontal linear solver. The linear solver is fully parallel, working with substructures determined by a parallel graph partitioner. The parallel performance of the nonlinear parallel algorithm is demonstrated using damage localization problems for two and three dimensional crack models. Convergence studies of the predicted local variation of damage at the crack tip are performed using discretizations with up to one million degrees of freedom. Implementation issues for damage localization problems are discussed, including the selection of independent variables and the use of Riks continuation method.We acknowledge the financial support of Ministry of Science and Technology by National Research Laboratory program. (grant number 00-N-NL-01-C-026)     

A monolithic geometric multigrid solver for fluid‐structure interactions in ALE formulation

We present a monolithic geometric multigrid solver for fluid‐structure interaction problems in Arbitrary Lagrangian Eulerian coordinates. The coupled dynamics of an incompressible fluid with nonlinear hyperelastic solids gives rise to very large and ill‐conditioned systems of algebraic equations. Direct solvers usually are out of question because of memory limitations, and standard coupled iterative solvers are seriously affected by the bad condition number of the system matrices. The use of partitioned preconditioners in Krylov subspace iterations is an option, but the convergence will be limited by the outer partitioning. Our proposed solver is based on a Newton linearization of the fully monolithic system of equations, discretized by a Galerkin finite element method. Approximation of the linearized systems is based on a monolithic generalized minimal residual method iteration, preconditioned by a geometric multigrid solver. The special character of fluid‐structure interactions is accounted for by a partitioned scheme within the multigrid smoother only. Here, fluid and solid field are segregated as Dirichlet–Neumann coupling. We demonstrate the efficiency of the multigrid iteration by analyzing 2d and 3d benchmark problems. While 2d problems are well manageable with available direct solvers, challenging 3d problems highly benefit from the resulting multigrid solver. Copyright © 2015 John Wiley & Sons, Ltd.     

A new parallel sparse direct solver: Presentation and numerical experiments in large‐scale structural mechanics parallel computing

The main purpose of this work is to present a new parallel direct solver: Dissection solver. It is based on LU factorization of the sparse matrix of the linear system and allows to detect automatically and handle properly the zero‐energy modes, which are important when dealing with DDM. A performance evaluation and comparisons with other direct solvers (MUMPS, DSCPACK) are also given for both sequential and parallel computations. Results of numerical experiments with a two‐level parallelization of large‐scale structural analysis problems are also presented: FETI is used for the global problem parallelization and Dissection for the local multithreading. In this framework, the largest problem we have solved is of an elastic solid composed of 400 subdomains running on 400 computation nodes (3200 cores) and containing about 165 millions dof. The computation of one single iteration consumes less than 20 min of CPU time. Several comparisons to MUMPS are given for the numerical computation of large‐scale linear systems on a massively parallel cluster: performances and weaknesses of this new solver are highlighted. Copyright © 2011 John Wiley & Sons, Ltd.     

An efficient out‐of‐core multifrontal solver for large‐scale unsymmetric element problems

In many applications where the efficient solution of large sparse linear systems of equations is required, a direct method is frequently the method of choice. Unfortunately, direct methods have a potentially severe limitation: as the problem size grows, the memory needed generally increases rapidly. However, the in‐core memory requirements can be limited by storing the matrix and its factors externally, allowing the solver to be used for very large problems. We have designed a new out‐of‐core package for the large sparse unsymmetric systems that arise from finite‐element problems. The code, which is called HSL _MA78 , implements a multifrontal algorithm and achieves efficiency through the use of specially designed code for handling the input/output operations and efficient dense linear algebra kernels. These kernels, which are available as a separate package called HSL _MA74 , use high‐level BLAS to perform the partial factorization of the frontal matrices and offer both threshold partial and rook pivoting. In this paper, we describe the design of HSL _MA78 and explain its user interface and the options it offers. We also describe the algorithms used by HSL _MA74 and illustrate the performance of our new codes using problems from a range of practical applications. Copyright © 2008 John Wiley & Sons, Ltd.     

SOLUTION OF LARGE LINEAR SYSTEMS ON PIPELINED SIMD MACHINES

We developed a direct out-of-core solver for dense non-symmetric linear systems of ‘arbitrary’ size N×N. The algorithm fully employs the Basic Linear Algebra Subprograms (BLAS), and can therefore easily be adapted to different computer architectures by using the corresponding optimized routines. We used blocked versions of left-looking and right-looking variants of LU decomposition to perform most of the operations in Level 3 BLAS, to reduce the number of I/O operations and to minimize the CPU time usage. The storage requirements of the algorithm are only 2N×NB data elements where NB≪N. Depending on the sustained floating point performance and the sustained I/O rate of the given hardware, we derived formulas that allow for choosing optimal values of NB to balance between CPU time and I/O time. We tested the algorithm by means of linear systems derived from 3D-BEM for strongly and weakly singular integral equations and from interpolation problems for scattered data on closed surfaces in ℝ³. It took only about 2⋅5 CPU minutes on a 5 GFLOPS vector computer SNI S600/20 to solve a linear system of size 10000, which corresponds to a performance of 4⋅3 GFLOPS; a value of NB=650 gives a reasonable I/O time and the necessary main storage size is about 13 Mwords. In addition, we compared the algorithm with (1) an out-of-core version of GMRES and (2) a wavelet transform followed by in-core GMRES after thresholding. At least for boundary integral equations of classical boundary value problems of potential theory, the out-of-core version of GMRES is superior to the direct out-of-core solver and the wavelet transform since the algorithm converged after at most 5 iteration steps. It took about 17 s to solve a system with 8192 unknowns compared with 146 s for direct out-of-core and 402 s for wavelet transform followed by in-core GMRES. © 1997 by John Wiley & Sons, Ltd.     

Two solution strategies to improve the computational performance of sequentially linear analysis for quasi-brittle structures

Sequentially linear analysis (SLA), an event-by-event procedure for finite element (FE) simulation of quasi-brittle materials, is based on sequentially identifying a critical integration point in the FE model, to reduce its strength and stiffness, and the corresponding critical load multiplier (λ_crit), to scale the linear analysis results. In this article, two strategies are proposed to efficiently reuse previous stiffness matrix factorisations and their corresponding solutions in subsequent linear analyses, since the global system of linear equations representing the FE model changes only locally. The first is based on a direct solution method in combination with the Woodbury matrix identity, to compute the inverse of a low-rank corrected stiffness matrix relatively cheaply. The second is a variation of the traditional incomplete LU preconditioned conjugate gradient method, wherein the preconditioner is the complete factorisation of a previous analysis step's stiffness matrix. For both the approaches, optimal points at which the factorisation is recomputed are determined such that the total analysis time is minimised. Comparison and validation against a traditional parallel direct sparse solver, with regard to a two-dimensional (2D) and three-dimensional (3D) benchmark study, illustrates the improved performance of the Woodbury-based direct solver over its counterparts, especially for large 3D problems.     

An algorithm for the matrix‐free solution of quasistatic frictional contact problems

A contact enforcement algorithm has been developed for matrix‐free quasistatic finite element techniques. Matrix‐free (iterative) solution algorithms such as non‐linear conjugate gradients (CG) and dynamic relaxation (DR) are desirable for large solid mechanics applications where direct linear equation solving is prohibitively expensive, but in contrast to more traditional Newton–Raphson and quasi‐Newton iteration strategies, the number of iterations required for convergence is typically of the same order as the number of degrees of freedom of the model. It is therefore crucial that each of these iterations be inexpensive to per‐form, which is of course the essence of a matrix free method. In applying such methods to contact problems we emphasize here two requirements: first, that the treatment of the contact should not make an average equilibrium iteration considerably more expensive; and second, that the contact constraints should be imposed in such a way that they do not introduce spurious energy that acts against the iterative solver. These practical concerns are utilized to develop an iterative technique for accurate constraint enforcement that is suitable for non‐linear conjugate gradient and dynamic relaxation iterative schemes. Copyright © 1999 John Wiley & Sons, Ltd.     

Fully implicit Lagrange–Newton–Krylov–Schwarz algorithms for boundary control of unsteady incompressible flows

We develop a parallel fully implicit domain decomposition algorithm for solving optimization problems constrained by time‐dependent nonlinear partial differential equations. In particular, we study the boundary control of unsteady incompressible Navier–Stokes equations. After an implicit discretization in time, a fully coupled sparse nonlinear optimization problem needs to be solved at each time step. The class of full space Lagrange–Newton–Krylov–Schwarz algorithms is used to solve the sequence of optimization problems. Among optimization algorithms, the fully implicit full space approach is considered to be the easiest to formulate and the hardest to solve. We show that Lagrange–Newton–Krylov–Schwarz, with a one‐level restricted additive Schwarz preconditioner, is an efficient class of methods for solving these hard problems. To demonstrate the scalability and robustness of the algorithm, we consider several problems with a wide range of Reynolds numbers and time step sizes, and we present numerical results for large‐scale calculations involving several million unknowns obtained on machines with more than 1000 processors. Copyright © 2012 John Wiley & Sons, Ltd.     

CMRH method as iterative solver for boundary element acoustic systems

Disceretization of boundary integral equations leads to complex and fully populated linear systems. One inherent drawback of the boundary element method (BEM) is that, the dense linear system has to be constructed and solved for each frequency. For large-scale problems, BEM can be more efficient by improving the construction and solution phases of the linear system. For these problems, the application of common direct solver is inefficient. In this paper, the corresponding linear systems are solved more efficiently than common direct solvers by using the iterative technique called CMRH (Changing Minimal Residual method based on Hessenberg process). In this method, the generation of the basis vectors of the Krylov subspace is based on the Hessenberg process instead of Arnoldi's one that the most known GMRES (Generalized Minimal RESidual) solver uses. Compared to GMRES, the storage requirements are considerably reduced in CMRH.     

Theoretical comparison of the FETI and algebraically partitioned FETI methods,and performance comparisons with a direct sparse solver

In this paper, we prove that the Algebraic A‐FETI method corresponds to one particular instance of the original one‐level FETI method. We also report on performance comparisons on an Origin 2000 between the one‐ and two‐level FETI methods and an optimized sparse solver, for two industrial applications: the stress analysis of a thin shell structure, and that of a three‐dimensional structure modelled by solid elements. These comparisons suggest that for topologically two‐dimensional problems, sparse solvers are effective when the number of processors is relatively small. They also suggest that for three‐dimensional applications, scalable domain decomposition methods such as FETI deliver a superior performance on both sequential and parallel hardware configurations. Copyright © 1999 John Wiley & Sons, Ltd.     

A modified Newton‐type Koiter‐Newton method for tracing the geometrically nonlinear response of structures

The Koiter‐Newton (KN) method is a combination of local multimode polynomial approximations inspired by Koiter's initial postbuckling theory and global corrections using the standard Newton method. In the original formulation, the local polynomial approximation, called a reduced‐order model, is used to make significantly more accurate predictions compared to the standard linear prediction used in conjunction with Newton method. The correction to the exact equilibrium path relied exclusively on Newton‐Raphson method using the full model. In this paper, we proposed a modified Newton‐type KN method to trace the geometrically nonlinear response of structures. The developed predictor‐corrector strategy is applied to each predicted solution of the reduced‐order model. The reduced‐order model can be used also in the correction phase, and the exact full nonlinear model is applied only to calculate force residuals. Remainder terms in both the displacement expansion and the reduced‐order model are well considered and constantly updated during correction. The same augmented finite element model system is used for both the construction of the reduced‐order model and the iterations for correction. Hence, the developed method can be seen as a particular modified Newton method with a constant iteration matrix over the single KN step. This significantly reduces the computational cost of the method. As a side product, the method has better error control, leading to more robust step size adaptation strategies. Numerical results demonstrate the effectiveness of the method in treating nonlinear buckling problems.     

Comparison of fully implicit and IMPES formulations for simulation of water injection in fractured and unfractured media

In this work, we compare the fully implicit (FI) and implicit pressure‐explicit saturation (IMPES) formulations for the simulation of water injection in fractured media. The system of partial differential equations is discretized within the discrete‐fracture framework using a control‐volume method. A unique feature of the methodology is that there is no need for the computation of matrix–fracture transfer terms. The non‐linear system of equations resulting from the FI formulation is solved with state‐of‐the‐art Newton and tensor methods. Direct and Krylov iterative methods are employed to solve the system resulting from the Newton linearization. The performance of the FI and IMPES formulations is compared with numerical testing. Results show that the contrast between matrix and fracture properties affects the performance of both IMPES and FI formulations and that the tensor method outperforms all the Newton solvers for the near‐singular Jacobian matrix. Copyright © 2006 John Wiley & Sons, Ltd.     

A Generalised Conjugate Residual method for the solution of non‐symmetric systems of equations with multiple right‐hand sides

This paper presents an iterative algorithm for solving non‐symmetric systems of equations with multiple right‐hand sides. The algorithm is an extension of the Generalised Conjugate Residual method (GCR) and combines the advantages of a direct solver with those of an iterative solver: it does not have to restart from scratch for every right‐hand side, it tends to require less memory than a direct solver, and it can be implemented efficiently on a parallel computer. We will show that the extended GCR algorithm can be competitive with a direct solver when running on a single processor. We will also show that the algorithm performs well on a Cray T3E parallel computer. Copyright © 1999 John Wiley & Sons, Ltd.     

An adaptive radial basis algorithm (ARBF) for expensive black-box mixed-integer constrained global optimization

Response surface methods based on kriging and radial basis function (RBF) interpolation have been successfully applied to solve expensive, i.e. computationally costly, global black-box nonconvex optimization problems. In this paper we describe extensions of these methods to handle linear, nonlinear, and integer constraints. In particular, algorithms for standard RBF and the new adaptive RBF (ARBF) are described. Note, however, while the objective function may be expensive, we assume that any nonlinear constraints are either inexpensive or are incorporated into the objective function via penalty terms. Test results are presented on standard test problems, both nonconvex problems with linear and nonlinear constraints, and mixed-integer nonlinear problems (MINLP). Solvers in the TOMLAB Optimization Environment () have been compared, specifically the three deterministic derivative-free solvers rbfSolve, ARBFMIP and EGO with three derivative-based mixed-integer nonlinear solvers, OQNLP, MINLPBB and MISQP, as well as the GENO solver implementing a stochastic genetic algorithm. Results show that the deterministic derivative-free methods compare well with the derivative-based ones, but the stochastic genetic algorithm solver is several orders of magnitude too slow for practical use. When the objective function for the test problems is costly to evaluate, the performance of the ARBF algorithm proves to be superior.     

Residual-based variational multiscale simulation of free surface flows

In this work we apply the residual-based variational multiscale method (RB-VMS) to the volume-of-fluid (VOF) formulation of free-surface flows. Using this technique we are able to solve such problems in a Large Eddy Simulation framework. This is a natural extension of our Navier–Stokes solver, which uses the RB-VMS finite element formulation, edge-based data structures, adaptive time step control, inexact Newton solvers and supports several parallel programming paradigms. The VOF interface capturing variable is advected using the computed coarse and fine scales velocity field. Thus, the RB-VMS technique can be readily applied to the free-surface solver with minor modifications on the implementation. We apply this technique to the solution of two problems where available data indicate complex free-surface behavior. Results are compared with numerical and experimental data and show that the present formulation can achieve good accuracy with minor impacts on computational efficiency.     

Analytic preconditioners for the electric field integral equation

Since the advent of the fast multipole method, large‐scale electromagnetic scattering problems based on the electric field integral equation (EFIE) formulation are generally solved by a Krylov iterative solver. A well‐known fact is that the dense complex non‐hermitian linear system associated to the EFIE becomes ill‐conditioned especially in the high‐frequency regime. As a consequence, this slows down the convergence rate of Krylov subspace iterative solvers. In this work, a new analytic preconditioner based on the combination of a finite element method with a local absorbing boundary condition is proposed to improve the convergence of the iterative solver for an open boundary. Some numerical tests precise the behaviour of the new preconditioner. Moreover, comparisons are performed with the analytic preconditioner based on the Calderòn's relations for integral equations for several kinds of scatterers. Copyright © 2004 John Wiley & Sons, Ltd.     

A fast boundary element based solver for localized inelastic deformations

We present a numerical method for the solution of nonlinear geomechanical problems involving localized deformation along shear bands and fractures. We leverage the boundary element method to solve for the quasi-static elastic deformation of the medium while rigid-plastic constitutive relations govern the behavior of displacement discontinuity (DD) segments capturing localized deformations. A fully implicit scheme is developed using a hierarchical approximation of the boundary element matrix. Combined with an adequate block preconditioner, this allows to tackle large problems via the use of an iterative solver for the solution of the tangent system. Several two-dimensional examples of the initiation and growth of shear-bands and tensile fractures illustrate the capabilities and accuracy of this technique. The method does not exhibit any mesh dependency associated with localization provided that (i) the softening length-scale is resolved and (ii) the plane of localized deformations is discretized a priori using DD segments.     

A fast implementation of the FETI‐DP method: FETI‐DP‐RBS‐LNA and applications on large scale problems with localized non‐linearities

As parallel and distributed computing gradually becomes the computing standard for large scale problems, the domain decomposition method (DD) has received growing attention since it provides a natural basis for splitting a large problem into many small problems, which can be submitted to individual computing nodes and processed in a parallel fashion. This approach not only provides a method to solve large scale problems that are not solvable on a single computer by using direct sparse solvers but also gives a flexible solution to deal with large scale problems with localized non‐linearities. When some parts of the structure are modified, only the corresponding subdomains and the interface equation that connects all the subdomains need to be recomputed. In this paper, the dual–primal finite element tearing and interconnecting method (FETI‐DP) is carefully investigated, and a reduced back‐substitution (RBS) algorithm is proposed to accelerate the time‐consuming preconditioned conjugate gradient (PCG) iterations involved in the interface problems. Linear–non‐linear analysis (LNA) is also adopted for large scale problems with localized non‐linearities based on subdomain linear–non‐linear identification criteria. This combined approach is named as the FETI‐DP‐RBS‐LNA algorithm and demonstrated on the mechanical analyses of a welding problem. Serial CPU costs of this algorithm are measured at each solution stage and compared with that from the IBM Watson direct sparse solver and the FETI‐DP method. The results demonstrate the effectiveness of the proposed computational approach for simulating welding problems, which is representative of a large class of three‐dimensional large scale problems with localized non‐linearities. Copyright © 2005 John Wiley & Sons, Ltd.