High-order Discretization Research Articles

Fully implicit higher-order schemes applied to polymer flooding

In water-based EOR methods, active chemical or biological substances are added to modify the physical properties of the fluids or/and the porous media at the interface between oil and water. The resulting displacement processes are governed by complex interplays between the transport of chemical substances, which is largely linear and highly affected by numerical diffusion, and how these substances affect the flow by changing the properties of the fluids and the surrounding rock. These effects are highly nonlinear and highly sensitive to threshold parameters that determine sharp transitions between regions of very different behavior. Unresolved simulation can therefore lead to misleading predictions of injectivity and recovery profiles. Use of higher-order spatial discretization schemes have been proposed by many authors as a means to reduce numerical diffusion and grid-orientation effects. Most higher-order simulators reported in the literature are based on explicit time stepping, and only a few are implicit. One reason that fully implicit formulations are not widely used might be that it becomes quite involved to compute the necessary linearizations for modern high-resolution discretizations of TVD and WENO type. Herein, we solve this problem by using automatic differentiation. We also demonstrate that using lagged evaluation of slope limiters and WENO weights alleviates the nonlinearity of the discrete systems and improves the computational efficiency, without having an adverse effect on the stability and accuracy of the higher-resolution schemes. As an example of EOR, we consider polymer flooding, which involves complex and adverse phenomena like adsorption in the rock, degradation and in-situ chemical reactions, shear thinning/thickening, dead pore space, etc. Using a few idealized test cases, we compare and contrast explicit and fully implicit time stepping for a variety of high and low-resolution spatial discretizations.

Application of a high-order CFD harmonic balance method to nonlinear aeroelasticity

An Aeroelastic-Harmonic Balance (A-HB) formulation of the Euler flow equations using a high-order spatial discretization scheme coupled with structural dynamic equations is proposed. The main objective of this new approach is to dramatically reduce the computational cost required to predict unsteady, periodic problems such as limit cycle oscillations (LCO). To this end, a new solver based on the Monotonicity Preserving limiter together with the AUSM+-up flux function is developed for the harmonic balance equations. The use of high-order CFD schemes allows the reduction of the number of degrees of freedom required to achieve a given desired accuracy, with respect to lower order schemes. In this paper, the reduction in degrees of freedom of the fluid system is exploited in the context of a CFD based Harmonic-Balance framework using a frequency updating procedure to determine the limit cycle conditions. The standard A-HB methodology has shown over one order of magnitude speed-up over time-marching methods; by employing the proposed high-order scheme in conjunction with coarser grids, the LCO computational time is halved without compromising accuracy.

Stochastic Galerkin approximation of the Reynolds equation with irregular film thickness

We consider the approximation of the Reynolds equation with an uncertain film thickness. The resulting stochastic partial differential equation is solved numerically by the stochastic Galerkin finite element method with high-order discretizations both in spatial and stochastic domains. We compute the pressure field of a journal bearing in various numerical examples that demonstrate the effectiveness and versatility of the approach. The results suggest that the stochastic Galerkin method is capable of supporting design when manufacturing imperfections are the main sources of uncertainty.

Uniform and high-order discretization schemes for Sturm\u2013Liouville problems via Fer streamers

The current paper concerns the uniform and high-order discretization of the novel approach to the computation of Sturm–Liouville problems via Fer streamers, put forth in Ramos and Iserles (Numer. Math. 131(3), 541—565 2015). In particular, the discretization schemes are shown to enjoy large step sizes uniform over the entire eigenvalue range and tight error estimates uniform for every eigenvalue. They are made explicit for global orders 4,7,10. In addition, the present paper provides total error estimates that quantify the interplay between the truncation and the discretization in the approach by Fer streamers.

Contact-aware simulations of particulate Stokesian suspensions

We present an efficient, accurate, and robust method for simulation of dense suspensions of deformable and rigid particles immersed in Stokesian fluid in two dimensions. We use a well-established boundary integral formulation for the problem as the foundation of our approach. This type of formulation, with a high-order spatial discretization and an implicit and adaptive time discretization, have been shown to be able to handle complex interactions between particles with high accuracy. Yet, for dense suspensions, very small time-steps or expensive implicit solves as well as a large number of discretization points are required to avoid non-physical contact and intersections between particles, leading to infinite forces and numerical instability.Our method maintains the accuracy of previous methods at a significantly lower cost for dense suspensions. The key idea is to ensure interference-free configuration by introducing explicit contact constraints into the system. While such constraints are unnecessary in the formulation, in the discrete form of the problem, they make it possible to eliminate catastrophic loss of accuracy by preventing contact explicitly.Introducing contact constraints results in a significant increase in stable time-step size for explicit time-stepping, and a reduction in the number of points adequate for stability.

A space-time adaptive method for reservoir flows: formulation and one-dimensional application

This paper presents a space-time adaptive framework for solving porous media flow problems, with specific application to reservoir simulation. A fully unstructured mesh discretization of space and time is used instead of a conventional time-marching approach. A space-time discontinuous Galerkin finite element method is employed to achieve a high-order discretization on the anisotropic, unstructured meshes. Anisotropic mesh adaptation is performed to reduce the error of a specified output of interest, by using a posteriori error estimates from the dual-weighted residual method to drive a metric-based mesh optimization algorithm. The space-time adaptive method is tested on a one-dimensional two-phase flow problem, and is found to be more efficient in terms of computational cost (degrees-of-freedom and total runtime) required to achieve a specified output error level, when compared to a conventional first-order time-marching finite volume method and the space-time discontinuous Galerkin method on structured meshes.

Implicit mesh discontinuous Galerkin methods and interfacial gauge methods for high-order accurate interface dynamics, with applications to surface tension dynamics, rigid body fluid–structure interaction, and free surface flow: Part II

In this two-part paper, a high-order accurate implicit mesh discontinuous Galerkin (dG) framework is developed for fluid interface dynamics, facilitating precise computation of interfacial fluid flow in evolving geometries. The framework uses implicitly defined meshes—wherein a reference quadtree or octree grid is combined with an implicit representation of evolving interfaces and moving domain boundaries—and allows physically prescribed interfacial jump conditions to be imposed or captured with high-order accuracy. Part one discusses the design of the framework, including: (i) high-order quadrature for implicitly defined elements and faces; (ii) high-order accurate discretisation of scalar and vector-valued elliptic partial differential equations with interfacial jumps in ellipticity coefficient, leading to optimal-order accuracy in the maximum norm and discrete linear systems that are symmetric positive (semi)definite; (iii) the design of incompressible fluid flow projection operators, which except for the influence of small penalty parameters, are discretely idempotent; and (iv) the design of geometric multigrid methods for elliptic interface problems on implicitly defined meshes and their use as preconditioners for the conjugate gradient method. Also discussed is a variety of aspects relating to moving interfaces, including: (v) dG discretisations of the level set method on implicitly defined meshes; (vi) transferring state between evolving implicit meshes; (vii) preserving mesh topology to accurately compute temporal derivatives; (viii) high-order accurate reinitialisation of level set functions; and (ix) the integration of adaptive mesh refinement.In part two, several applications of the implicit mesh dG framework in two and three dimensions are presented, including examples of single phase flow in nontrivial geometry, surface tension-driven two phase flow with phase-dependent fluid density and viscosity, rigid body fluid–structure interaction, and free surface flow. A class of techniques known as interfacial gauge methods is adopted to solve the corresponding incompressible Navier–Stokes equations, which, compared to archetypical projection methods, have a weaker coupling between fluid velocity, pressure, and interface position, and allow high-order accurate numerical methods to be developed more easily. Convergence analyses conducted throughout the work demonstrate high-order accuracy in the maximum norm for all of the applications considered; for example, fourth-order spatial accuracy in fluid velocity, pressure, and interface location is demonstrated for surface tension-driven two phase flow in 2D and 3D. Specific application examples include: vortex shedding in nontrivial geometry, capillary wave dynamics revealing fine-scale flow features, falling rigid bodies tumbling in unsteady flow, and free surface flow over a submersed obstacle, as well as high Reynolds number soap bubble oscillation dynamics and vortex shedding induced by a type of Plateau–Rayleigh instability in water ripple free surface flow. These last two examples compare numerical results with experimental data and serve as an additional means of validation; they also reveal physical phenomena not visible in the experiments, highlight how small-scale interfacial features develop and affect macroscopic dynamics, and demonstrate the wide range of spatial scales often at play in interfacial fluid flow.

Implicit mesh discontinuous Galerkin methods and interfacial gauge methods for high-order accurate interface dynamics, with applications to surface tension dynamics, rigid body fluid–structure interaction, and free surface flow: Part I

In this two-part paper, a high-order accurate implicit mesh discontinuous Galerkin (dG) framework is developed for fluid interface dynamics, facilitating precise computation of interfacial fluid flow in evolving geometries. The framework uses implicitly defined meshes—wherein a reference quadtree or octree grid is combined with an implicit representation of evolving interfaces and moving domain boundaries—and allows physically prescribed interfacial jump conditions to be imposed or captured with high-order accuracy. Part one discusses the design of the framework, including: (i) high-order quadrature for implicitly defined elements and faces; (ii) high-order accurate discretisation of scalar and vector-valued elliptic partial differential equations with interfacial jumps in ellipticity coefficient, leading to optimal-order accuracy in the maximum norm and discrete linear systems that are symmetric positive (semi)definite; (iii) the design of incompressible fluid flow projection operators, which except for the influence of small penalty parameters, are discretely idempotent; and (iv) the design of geometric multigrid methods for elliptic interface problems on implicitly defined meshes and their use as preconditioners for the conjugate gradient method. Also discussed is a variety of aspects relating to moving interfaces, including: (v) dG discretisations of the level set method on implicitly defined meshes; (vi) transferring state between evolving implicit meshes; (vii) preserving mesh topology to accurately compute temporal derivatives; (viii) high-order accurate reinitialisation of level set functions; and (ix) the integration of adaptive mesh refinement.In part two, several applications of the implicit mesh dG framework in two and three dimensions are presented, including examples of single phase flow in nontrivial geometry, surface tension-driven two phase flow with phase-dependent fluid density and viscosity, rigid body fluid–structure interaction, and free surface flow. A class of techniques known as interfacial gauge methods is adopted to solve the corresponding incompressible Navier–Stokes equations, which, compared to archetypical projection methods, have a weaker coupling between fluid velocity, pressure, and interface position, and allow high-order accurate numerical methods to be developed more easily. Convergence analyses conducted throughout the work demonstrate high-order accuracy in the maximum norm for all of the applications considered; for example, fourth-order spatial accuracy in fluid velocity, pressure, and interface location is demonstrated for surface tension-driven two phase flow in 2D and 3D. Specific application examples include: vortex shedding in nontrivial geometry, capillary wave dynamics revealing fine-scale flow features, falling rigid bodies tumbling in unsteady flow, and free surface flow over a submersed obstacle, as well as high Reynolds number soap bubble oscillation dynamics and vortex shedding induced by a type of Plateau–Rayleigh instability in water ripple free surface flow. These last two examples compare numerical results with experimental data and serve as an additional means of validation; they also reveal physical phenomena not visible in the experiments, highlight how small-scale interfacial features develop and affect macroscopic dynamics, and demonstrate the wide range of spatial scales often at play in interfacial fluid flow.

Numerical simulation of 1-D oil and water displacements in petroleum reservoirs using the spectral finite volume method

The study and development of high-order discretization techniques for the modeling and simulation of multiphase flows in petroleum reservoirs are still a challenge from the computational viewpoint due to the difficulty posed by some physical features such as heterogeneity and anisotropy of the medium, that are of paramount importance in this class of applications. In this paper, we present a preliminary investigation of the very high order Spectral Finite Volume Method (SFVM) for the discretization of the saturation equation that describes the 1-D Oil-Water displacement through porous media. In this equation, the hyperbolic term is discretized using a high order SFVM, the capillary term is approximated using a second order CVFD (Control Volume Finite Difference) method and the integration in time is carried out using a third-order Runge–Kutta method. On the other hand, to prevent numerical oscillations near shocks, that are typical from higher order schemes, a slope limiter procedure is used in the reconstruction stage. To evaluate the performance of the SFVM some simple model problems, including capillarity and gravity, are solved, yielding accurate and efficient results. The use of very high order formulations can be of utmost importance, particularly for multidimensional and real large-scale applications with complex wave structures.

Simulation of coupled structural-acoustic response with dynamic damage evolution

A novel time-domain method for simulating dynamic damage evolution in a coupled structural-acoustic system is presented. The system is derived via the theory of continuum damage mechanics, and incorporates standard damage evolution models, but is readily extendible to more exotic formulations. The overall solution method is staggered, solving for the dynamic damage evolution first with an explicit step, and then using the new values in the coupled computation of the structural-acoustic system. The spatial domain is discretized using a mixed finite element method, and the temporal domain is discretized with a higher-order implicit time discretization scheme. Efforts toward fully coupled verification of the solution algorithm are presented, as are validation studies for cases without evolving damage. Applications with evolving damage are presented, and present a novel first principles study of changes in the structural acoustic response to dynamically evolving damage in the structure. Special attention is given to brittle fracture. Examples of downstream usage of the evolving structural response are discussed in the concluding remarks.

Cahn–Hilliard on surfaces: A numerical study

The Cahn-Hilliard system has been used to describe a wide number of phase separation processes, from co-polymer systems to lipid membranes. In this work the convergence properties of a closest-point based scheme is investigated. In place of solving the original fourth-order system directly, two coupled second-order systems are solved. The system is solved using an incomplete Schur-decomposition as a preconditioner. The results indicate that with a sufficiently high-order time discretization the method only depends on the underlying spatial resolution.

Ws,p-approximation properties of elliptic projectors on polynomial spaces, with application to the error analysis of a Hybrid High-Order discretisation of Leray–Lions problems

In this work, we prove optimal [Formula: see text]-approximation estimates (with [Formula: see text]) for elliptic projectors on local polynomial spaces. The proof hinges on the classical Dupont–Scott approximation theory together with two novel abstract lemmas: An approximation result for bounded projectors, and an [Formula: see text]-boundedness result for [Formula: see text]-orthogonal projectors on polynomial subspaces. The [Formula: see text]-approximation results have general applicability to (standard or polytopal) numerical methods based on local polynomial spaces. As an illustration, we use these [Formula: see text]-estimates to derive novel error estimates for a Hybrid High-Order (HHO) discretisation of Leray–Lions elliptic problems whose weak formulation is classically set in [Formula: see text] for some [Formula: see text]. This kind of problems appears, e.g. in the modelling of glacier motion, of incompressible turbulent flows, and in airfoil design. Denoting by [Formula: see text] the meshsize, we prove that the approximation error measured in a [Formula: see text]-like discrete norm scales as [Formula: see text] when [Formula: see text] and as [Formula: see text] when [Formula: see text].

On the accuracy and efficiency of point-value multiresolution algorithms for solving scalar wave and Euler equations

This paper concerns the use of a point-value multiresolution algorithm and its extension to three-dimensional hyperbolic conservation laws. The proposed method is applied to a high-order finite-differences discretization with an explicit time integration. The fluxes are evaluated on the adaptive grid using a fifth-order high-resolution shock capturing scheme based on a WENO solver, while the time is advanced using a third-order Runge–Kutta scheme. The multiresolution prediction operators are presented for one-, two- and three-dimensional problems. To assess the efficiency and the accuracy of the method, a new tolerance-scale diagram is introduced. This diagram enables to properly choose the adequate value of the tolerance in order to maintain an optimal multiresolution quality. Numerical examples based on advection and Euler equations are carried out to show that the proposed method yields accurate results.

Variational boundary conditions based on the Nitsche method for fitted and unfitted isogeometric discretizations of the mechanically coupled Cahn–Hilliard equation

The primal variational formulation of the fourth-order Cahn–Hilliard equation requires C1-continuous finite element discretizations, e.g., in the context of isogeometric analysis. In this paper, we explore the variational imposition of essential boundary conditions that arise from the thermodynamic derivation of the Cahn–Hilliard equation in primal variables. Our formulation is based on the symmetric variant of Nitsche's method, does not introduce additional degrees of freedom and is shown to be variationally consistent. In contrast to strong enforcement, the new boundary condition formulation can be naturally applied to any mapped isogeometric parametrization of any polynomial degree. In addition, it preserves full accuracy, including higher-order rates of convergence, which we illustrate for boundary-fitted discretizations of several benchmark tests in one, two and three dimensions. Unfitted Cartesian B-spline meshes constitute an effective alternative to boundary-fitted isogeometric parametrizations for constructing C1-continuous discretizations, in particular for complex geometries. We combine our variational boundary condition formulation with unfitted Cartesian B-spline meshes and the finite cell method to simulate chemical phase segregation in a composite electrode. This example, involving coupling of chemical fields with mechanical stresses on complex domains and coupling of different materials across complex interfaces, demonstrates the flexibility of variational boundary conditions in the context of higher-order unfitted isogeometric discretizations.

Anderson Transition of Cold Atoms with Synthetic Spin-Orbit Coupling in Two-Dimensional Speckle Potentials.

We investigate the metal-insulator transition occurring in two-dimensional (2D) systems of noninteracting atoms in the presence of artificial spin-orbit interactions and a spatially correlated disorder generated by laser speckles. Based on a high order discretization scheme, we calculate the precise position of the mobility edge and verify that the transition belongs to the symplectic universality class. We show that the mobility edge depends strongly on the mixing angle between Rashba and Dresselhaus spin-orbit couplings. For equal couplings a non-power-law divergence is found, signaling the crossing to the orthogonal class, where such a 2D transition is forbidden.

SPARC: Accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory: Extended systems

As the second component of SPARC (Simulation Package for Ab-initio Real-space Calculations), we present an accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory (DFT) for extended systems. Specifically, employing a local formulation of the electrostatics, the Chebyshev polynomial filtered self-consistent field iteration, and a reformulation of the non-local force component, we develop a finite-difference framework wherein both the energy and atomic forces can be efficiently calculated to within desired accuracies in DFT. We demonstrate using a wide variety of materials systems that SPARC achieves high convergence rates in energy and forces with respect to spatial discretization to reference plane-wave result; exponential convergence in energies and forces with respect to vacuum size for slabs and wires; energies and forces that are consistent and display negligible ‘egg-box’ effect; accurate properties of crystals, slabs, and wires; and negligible drift in molecular dynamics simulations. We also demonstrate that the weak and strong scaling behavior of SPARC is similar to well-established and optimized plane-wave implementations for systems consisting up to thousands of electrons, but with a significantly reduced prefactor. Overall, SPARC represents an attractive alternative to plane-wave codes for performing DFT simulations of extended systems. Program summaryProgram title: SPARCCatalogue identifier: AFBR_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AFBR_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: GNU GPL v3No. of lines in distributed program, including test data, etc.: 93822No. of bytes in distributed program, including test data, etc.: 1386659Distribution format: tar.gzProgramming language: C/C++.Computer: Any system with C/C++ compiler.Operating system: Linux.RAM: Problem dependent. Ranges from 80 GB to 800 GB for a system with 2500 electrons.Classification: 7.3.External routines: PETSc 3.5.3 (http://www.mcs.anl.gov/petsc), MKL 11.2 (https://software.intel.com/en-us/intel-mkl), and MVAPICH2 2.1 (http://mvapich.cse.ohio-state.edu/).Does the new version supersede the previous version?: YesNature of problem: Calculation of the static and dynamic properties of isolated and extended systems in the framework of Kohn–Sham Density Functional Theory (DFT).Solution method:High-order finite-difference discretization. Local reformulation of the electrostatics in terms of the electrostatic potential and pseudocharge densities. Application of Bloch-periodic and zero-Dirichlet boundary conditions on the orbitals in the direction of periodicity and vacuum, respectively. Application of periodic and Dirichlet boundary conditions on the electrostatic potential in the direction of periodicity and vacuum, respectively. Integration over the Brillouin zone for extended systems using the Monkhorst–Pack grid. Calculation of the electronic ground-state using the Chebyshev polynomial filtered self-consistent field iteration in conjunction with Anderson based extrapolation/mixing schemes. Reformulation of the non-local component of the force. Geometry optimization using the Polak–Ribiere variant of non-linear conjugate gradients with secant line search. NVE molecular dynamics using the leapfrog method. Parallelization via domain decomposition and over Brillouin zone integration.Reasons for new version:To enable the study of extended systems like crystals, slabs, and wires using SPARC.Summary of revisions:Incorporated the ability to study the static and dynamic properties of crystals, slabs, and wires.Restrictions:System size less than ∼4000 electrons. Local Density Approximation (LDA). Troullier–Martins pseudopotentials without relativistic or non-linear core corrections. Domain has to be cuboidal.Running time:Problem dependent. Timing results for selected examples provided in the paper.

On the arithmetic intensity of high-order finite-volume discretizations for hyperbolic systems of conservation laws

It has been conjectured that higher-order discretizations for partial differential equations will have advantages over the lower-order counterparts commonly used today. The reasoning is that the increase in arithmetic operations will be more than offset by the reduction in data transfers and the increase in concurrent floating-point units. To evaluate this conjecture, the arithmetic intensity of a class of high-order finite-volume discretizations for hyperbolic systems of conservation laws is theoretically analyzed for spatial discretizations from orders three through eight in arbitrary dimensions. Three cache models are considered: the limiting cases of no cache and infinite cache as well as a finite-sized cache model. Models are validated experimentally by measuring floating-point operations and data transfers on an IBM Blue Gene/Q node. Theory and experiments demonstrate that high-order finite-volume methods will be able to provide increases in arithmetic intensity that will be necessary to make better utilization of on-node floating-point capability.

Multigrid cell-centered techniques for high-order incompressible flow numerical solutions

A multigrid pressure correction scheme suitable for high order discretizations of the incompressible Navier–Stokes equations is developed and demonstrated. The pressure correction equation is discretized with fourth-order compact finite-difference approximations. Iterative methods based on multigrid techniques accelerate the most demanding part of the overall solution algorithm, which is the numerical solution of the arised large and sparse linear system. Geometrical multigrid methods, using partial semicoarsenig strategy and zebra line Gauss–Seidel relaxation, are employed to efficiently approximate the solution of the resulting algebraic linear system. Effects of various multigrid components on the pressure correction procedure are evaluated and new high-order transfer operators are developed for the case of cell-centered grids. Their convergence rates are also compared with commonly used intergrid transfer operators. Furthermore, numerically comparisons between different multigrid cycle approaches, such as V-, W- and F-cycle, are presented. The performance tests demonstrate that the new pressure correction approach significantly reduces the computational effort compared to single-grid algorithms. Furthermore, it is shown that the overall high order accuracy of the numerical method is retained in space and time with increasing Reynolds number.

Development of a hybrid parallel MCV-based high-order global shallow-water model

Utilization of high-order spatial discretizations is an important trend in developing global atmospheric models. As a competitive choice, the multi-moment constrained volume (MCV) method can achieve high accuracy while maintaining similar parallel scalability to classical finite volume methods. In this work, we introduce the development of a hybrid parallel MCV-based global shallow-water model on the cubed-sphere grid. Based on a sequential code, we perform parallelization on both the process and the thread levels. To enable process-level parallelism, we first decompose the six patches of the cubed-sphere in a same 2-D partition and then employ a conflict-free pipe-flow communication scheme for overlapping the halo exchange with computations. To further exploit the heterogeneous computing capacity of an Intel Xeon Phi accelerated supercomputer, we propose a guided panel-based inner---outer partition to distribute workload among the CPUs and the coprocessors. In addition to the above, thread-level parallelism along with various optimizations is done on both the multi-core CPU and the many-core accelerator. Numerical experiments are carried out to validate the correctness of the optimized parallel code and examine its parallel performance. Test results show that both the CPU-only and the hybrid codes scale well to hundreds of processes in terms of both the strong and weak scaling. In particular, the hybrid code can achieve a speedup of $$2.56\times $$2.56× as compared to the CPU-only version. In the largest run on a $$9216\,\times \,9216\,\times \,6$$9216×9216×6 mesh (1.5 billion unknowns), the hybrid code sustains an aggregative performance of 26.5 Tflops with 486 processes (33,534 cores).

Algebraic structure of vector fields in financial diffusion models and its applications

High-order discretization schemes of SDEs using free Lie algebra-valued random variables are introduced by Kusuoka [Adv. Math. Econ., 2004, 5, 69–83], [Adv. Math. Econ., 2013, 17, 71–120], Lyons–Victoir [Proc. R. Soc. Lond. Ser. A Math. Phys. Sci., 2004, 460, 169–198], Ninomiya–Victoir [Appl. Math. Finance, 2008, 15, 107–121] and Ninomiya–Ninomiya [Finance Stochast., 2009, 13, 415–443]. These schemes are called KLNV methods. They involve solving the flows of vector fields associated with SDEs and it is usually done by numerical methods. The authors have found a special Lie algebraic structure on the vector fields in the major financial diffusion models. Using this structure, we can solve the flows associated with vector fields analytically and efficiently. Numerical examples show that our method reduces the computation time drastically.