CPU Time Savings Research Articles

Optimal Combination of the Splitting–Linearizing Method to SSOR and SAOR for Solving the System of Nonlinear Equations

The symmetric successive overrelaxation (SSOR) and symmetric accelerated overrelaxation (SAOR) are conventional iterative methods for solving linear equations. In this paper, novel approaches are presented by combining a splitting–linearizing method with SSOR and SAOR for solving a system of nonlinear equations. The nonlinear terms are decomposed at two sides through a splitting parameter, which are linearized around the values at the previous step, obtaining a linear equation system at each iteration step. The optimal values of parameters are determined to minimize the reciprocal of the maximal projection, which are sought in preferred ranges using the golden section search algorithm. Numerical tests assess the performance of the developed methods, namely, the optimal splitting symmetric successive over-relaxation (OSSSOR), and the optimal splitting symmetric accelerated over-relaxation (OSSAOR). The chief advantages of the proposed methods are that they do not need to compute the inverse matrix at each iteration step, and the computed orders of convergence by OSSSOR and OSSAOR are between 1.5 and 5.61; they, without needing the inner iterations loop, converge very fast with saving CPU time to find the true solution with a high accuracy.

Model order reduction based on low-rank approximation for parameterized eigenvalue problems in structural dynamics

This study introduces a novel reduced-order modeling approach for parameterized eigenvalue problems. The proposed framework represents eigenvectors as a low-rank approximation and is implemented using an offline-online strategy. In the offline stage, nonlinear equations are formulated by imposing a stationary condition on the weighted Rayleigh quotient over sampled points and eigenvalues of interest. A variant of the Krylov subspace approach is applied to progressively identify the dominant response characteristics at sampled points. In the online stage, a reduced-order model is obtained through Galerkin projection, and eigensolutions at unsampled points are determined. The performance of the proposed framework is validated through numerical examples, where perturbation, least-angle regression, and Gaussian process are considered for comparison. The obtained results demonstrate that the proposed framework can accurately estimate eigensolutions, even when eigenvalues are adjacent to each other. In terms of efficiency, the a priori construction of the low-dimensional basis is possible without solving the full-order model, resulting in substantial CPU time savings during the offline computation process. Lastly, the applicability of the proposed framework is investigated through an uncertainty quantification example using numerical and experimental results.

Proper generalized decomposition-based iterative enrichment process combined with shooting method for steady-state forced response analysis of nonlinear dynamical systems

This paper presents a novel framework combining proper generalized decomposition (PGD) with the shooting method to determine the steady-state response of nonlinear dynamical systems upon a general periodic input. The proposed PGD approximates the response as a low-rank separated representation of the spatial and temporal dimensions. The Galerkin projection is employed to formulate the subproblem for each dimension, then the fixed-point iteration is applied. The subproblem for the spatial vector can be regarded as computing a set of reduced-order basis vectors, and the shooting problem projected onto the subspace spanned by these basis vectors is defined to obtain the temporal coefficients. From this procedure, the proposed framework replaces the complex nonlinear time integration of the full-order model with the series of solving simple iterative subproblems. The proposed framework is validated through two descriptive numerical examples considering the conventional linear normal mode method for comparison. The results show that the proposed shooting method based on PGD can accurately capture nonlinear characteristics within 10 modes, whereas linear modes cannot easily approximate these behaviors. In terms of computational efficiency, the proposed method enables CPU time savings of about one order of magnitude compared with the conventional shooting methods.

Investigation of Modified EMST Micromixing Model Performance on Lagrangian PDF Transported in Lifted Hydrogen/Air Fames

This article presents a numerical simulation of high hydrogen/air flames using the Lagrangian transported PDF method. This method enables the calculation of fluid composition changes resulting from convection and reaction without the need for modelling, while requiring modelling for molecular mixing. Consequently, the accuracy of calculations in this L-PDF method heavily relies on an accurate representation of the mixture model term. The Euclidean model, which provides a better description of physical mixing processes, is well-suited for modelling the molecular mixing term EMST. Additionally, the accuracy of this model depends on the value of the mixing constant, representing the ratio between the mechanical time scale and the scalar time scale. Two algebraic models for the mixing constant have been implemented in the computational code, employing a well-defined function to calculate this ratio for each cell. These models contribute to memory and CPU time savings. To account for turbulence and its interaction with physical phenomena, the RSM model is employed due to its ability to identify different areas of turbulent stresses. Hence, the primary objective of this study is to evaluate the capabilities of these algebraic models in predicting scalar fields within such flames. Overall, the predictions align well with experimental data, affirming the validity of these models.

Efficient sampling techniques for surrogate-based optimization with thermoelastic application

In many real-life models, the objective function may be very expensive to calculate or sometimes a black-box function. Therefore, a surrogate model is required to reduce the cost of the calculations. To build a robust surrogate model, a good sampling technique is needed. This paper suggests two new sampling techniques for surrogate-based optimization, depending on four random walk steps, which in turn give a good chance for building an accurate approximation for the original function. In addition to that, one of these new sampling techniques has a new feature that was not used before. The main advantage of our two suggested sampling techniques is that there is no optimization subproblem will be solved for the next sample point, which in turn saves CPU time. Furthermore, we applied our two sampling techniques to a 1D thermal shock problem for a half-space in the context of the fractional order theory of thermoelasticity to reduce the number of function evaluations and consequently CPU time. Laplace transform techniques are used to obtain the solution of the system in the frequency domain. Applying our algorithms to the few points obtained by using Laplace’s inverse, we can reach the optimal solution for all functions filed. The predictions of our techniques are discussed and compared with the classical method. 2D and 3D graphs are plotted for our numerical results. The temperature, displacement, and stress distributions are computed and represented graphically for different times and several fractional parameters.

Decomposition heuristics for multiobjective problems. The Food bank network redesign case

Large MO-MILP problems are often very hard to solve by exact methods. In this study we consider the Food bank network redesign (FBNR) problem to introduce three decompose-and-fix heuristics that successfully solve large real-life based instances of the problem. The FBNR problem is modeled as a multi-period, multi-product supply chain redesign problem that accounts for economic, environmental and social objectives in three distinct functions. Each new heuristic decomposes the FBNR problem into two MO-MILP problems, which are sequentially solved following the lexicographic concept. Decomposition observe the nature and terms of the decisions involved. The first MO-MILP problem concerns only the longer term decisions. After the corresponding decisions are fixed, the second MO-MILP problem referring to shorter term decisions is solved respecting the ranking of objectives followed for each solution obtained in the first problem. The difference among the heuristics lies in what is considered as longer or shorter term decisions. Besides solving large instances of the problem, which the exact method used could not do, CPU time savings generated by the heuristics range from 80% to 97% of the time required by the exact method. Compromise between the computational effort and the quality of solutions obtained by each heuristic is discussed. The solving methodology proposed in this study can be adapted to other large multiobjective problems.

Accelerating inverse inference of ensemble Kalman filter via reduced-order model trained using adaptive sparse observations

The implementation of data assimilation often struggles with low computational efficiency due to the complexity of high-fidelity numerical modeling and the large size of the observation vectors. To address this challenge, it is necessary to employ reduced-order strategies for both the full-order dynamical model and observation operator. Proper orthogonal decomposition (POD) based reduced-order model (ROM) and discrete empirical interpolation method (DEIM) based sparse observation identification are often considered as primary solutions. However, their effective combination requires the determination of high-quality POD reduced-order basis vectors. The current research focuses on incorporating the ROM and DEIM into a Kalman filter framework, resulting in a reduced-order Kalman filter with sparse observations. The key aspect is identifying the optimal POD basis, which can be achieved iteratively by introducing a relative root-mean-square-error (RMSE) for a near-optimal exploration of both the desired low-order POD basis and observation locations. The singular value decomposition (SVD) of the snapshots is performed in a low-dimensional subspace, leading to significant computational savings. This improved POD basis not only enhances the reconstruction of the ROM, but also helps build the DEIM observation matrix, resulting in simultaneous dimensionality reduction of the state and observation space. The effectiveness of this algorithm is demonstrated through the retrieval of initial conditions for one-dimensional tsunami wave equations with real data scenarios. Experiments are conducted for both smooth and nonsmoothed initial conditions. The results show that when the observation data is available at 3 and 6 spatial observation locations, our method only requires fewer vectors with the first 45 and 90 dominant modes derived using the RMSE, compared to as many as 210 and 252 derived from the relative energy criterion. Our method enables a low-modal ROM to effectively reconstruct the wave height field, leading to considerable prediction capability for potential long-wave dynamics with comparable accuracy. The CPU time savings achieved through our method is at least one to two orders of magnitude compared to the CPU time required by the full-order model.

A Robust Framework for Generating Adsorption Isotherms to Screen Materials for Carbon Capture.

To rank the performance of materials for a given carbon capture process, we rely on pure component isotherms from which we predict the mixture isotherms. For screening a large number of materials, we also increasingly rely on isotherms predicted from molecular simulations. In particular, for such screening studies, it is important that the procedures to generate the data are accurate, reliable, and robust. In this work, we develop an efficient and automated workflow for a meticulous sampling of pure component isotherms. The workflow was tested on a set of metal-organic frameworks (MOFs) and proved to be reliable given different guest molecules. We show that the coupling of our workflow with the Clausius-Clapeyron relation saves CPU time, yet enables us to accurately predict pure component isotherms at the temperatures of interest, starting from a reference isotherm at a given temperature. We also show that one can accurately predict the CO2 and N2 mixture isotherms using ideal adsorbed solution theory (IAST). In particular, we show that IAST is a more reliable numerical tool to predict binary adsorption uptakes for a range of pressures, temperatures, and compositions, as it does not rely on the fitting of experimental data, which typically needs to be done with analytical models such as dual-site Langmuir (DSL). This makes IAST a more suitable and general technique to bridge the gap between adsorption (raw) data and process modeling. To demonstrate this point, we show that the ranking of materials, for a standard three-step temperature swing adsorption (TSA) process, can be significantly different depending on the thermodynamic method used to predict binary adsorption data. We show that, for the design of processes that capture CO2 from low concentration (0.4%) streams, the commonly used methodology to predict mixture isotherms incorrectly assigns up to 33% of the materials as top-performing.

A Population-Based Iterated Greedy Algorithm for Distributed Assembly No-Wait Flow-Shop Scheduling Problem

This paper investigates a distributed assembly no-wait flow-shop scheduling problem (DANWFSP), which has important applications in manufacturing systems. The objective is to minimize the total flowtime. A mixed-integer linear programming model of DANWFSP with total flowtime criterion is proposed. A population-based iterated greedy algorithm (PBIGA) is presented to address the problem. A new constructive heuristic is presented to generate an initial population with high quality. For DANWFSP, an accelerated NR3 algorithm is proposed to assign jobs to the factories, which improves the efficiency of the algorithm and saves CPU time. To enhance the effectiveness of the PBIGA, the local search method and the destruction-construction mechanisms are designed for the product sequence and job sequence, respectively. A selection mechanism is presented to determine which individuals execute the local search method. An acceptance criterion is proposed to determine whether the offspring are adopted by the population. Last, the PBIGA and seven state-of-the-art algorithms are tested on 810 large-scale benchmark instances. The experimental results show that the presented PBIGA is an effective algorithm to address the problem and performs better than recently state-of-the-art algorithms compared in this paper.

Efficient quasi-brittle fracture simulations of concrete at mesoscale using micro CT images and a localizing gradient damage model

In this work, a localizing gradient damage model (LGDM) based on the generalized micromorphic theory is adopted to investigate the quasi-brittle fracture behaviour of concrete at mesoscale for the first time. Micro Computed Tomography (CT) image-based realistic concrete models are generated including aggregates, mortar, pores, and aggregate-mortar interfacial transition zones (ITZs). In the LGDM, an additional micro-force equilibrium equation is formulated with a damage-dependent interaction function to prevent spurious energetic interactions during the softening stage. The damage driving forces are evaluated as history variables through the maximum nonlocal equivalent strains at the integration points. The stiffness matrices and residual vectors for displacement and nonlocal equivalent strain fields are formulated through the conventional shape function and strain matrices of finite elements, by using the user-defined element (UEL) subroutine implemented with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm to improve numerical performance. Typical nonlinear fracture benchmarks are presented to investigate the variability of mesoscale fracture evolution, crack trajectories, and load–displacement curves as well as numerical performance. It is found that the LGDM significantly improves the computational efficiency over the conventional phase-field cohesive zone model (PFCZM) with 60% savings of CPU time, and holds great potential for multiscale fracture evaluation of quasi-brittle composite materials, through the high-fidelity representation of real materials’ microstructures, and flexible simulation of the complicated nonlinear fracture without the need for re-meshing or mesh enrichments.

An efficient model for predicting heat extraction from a multiple-well enhanced geothermal system in the presence of areal flow

To maximize the heat production from an enhanced geothermal system (EGS), multiple-well layouts have been proposed. In this paper, a semi-analytical model is developed for predicting heat extraction from a multiple-well system including injection and production wells, which are connected by a single planar aquifer with a constant thickness. By utilizing velocity potentials and streamline functions, analytical solutions for the temperature along streamlines in the Laplace space are obtained, followed by their numerical inversion to time-domain results. The model is verified through comparisons with existing results and those obtained by commercial package. Besides, the model is found to be very efficient in terms of seconds per run, compared to pure numerical models, in saving CPU time and providing acceptable results. After that, the thermal performance of eight different well layouts is compared to find the optimal well layout in terms of heat extraction for the EGS at Qiabuqia geothermal area in the Gonghe Basin, China. It is found that the magnitude and direction of the areal flow, and the well layout significantly influence the fluid flow and thermal behaviors. When the areal flow velocity is low (i.e. 10-5 m/s), the layout with 4 injection wells (IWs) and 5 production wells (PWs) exhibit the overall best output performance in terms of output thermal power (OTP) and total heat extracted (THE). However, if the areal flow velocity is high (i.e. 5 × 10-5 m/s), a dipole-well system, can maximize the OTP and THE if the two wells are placed parallel to the areal flow with the downstream IW and upstream PW.

Adaptive Multilevel Nonuniform Grid Algorithm for the Accelerated Analysis of Composite Metallic–Dielectric Radomes

An adaptive multilevel nonuniform grid (MLNG) algorithm is developed for the accelerated computation of fields radiated through composite metallic–dielectric radomes as well as antenna-radome interactions. The MLNG approach is applied to the mixed potential formulation of the coupled surface and volume electric field integral equations. The radome is decomposed into a hierarchy of subdomains (SDs) by an adaptive algorithm that closely follows the radome geometry, allowing significant savings in memory and CPU time. In the MLNG algorithm, only local generalized impedance matrices of the finest level SDs are evaluated. Far-zone potentials and fields are indirectly evaluated through a multilevel aggregation involving phase-and amplitude-compensated interpolation on nonuniform grids (NGs), requiring considerably fewer calculations as compared with the classical method of moments (MoM). The MLNG algorithm is incorporated in a preconditioned iterative solver. Accuracy as well as memory consumption and computation times of the algorithm are studied on realistic examples. The radome effect on the antenna input impedance and electric current density distribution is demonstrated. The method is validated by comparison with a commercial MoM software and shown to exhibit a computational complexity (CC) of O(NlogN), N being the number of unknowns.

A general and efficient numerical solution of reactive transport with multirate mass transfer

The presence of low permeability regions within porous media impacts solute transport and the distribution of species concentrations. Therefore, (bio)chemical reactions are equally affected. Multirate Mass Transfer (MRMT) models can be used to represent this anomalous transport process. MRMT conceptualizes the medium as a set of multiple continua: one mobile zone and multiple immobile zones. It simulates species transport in mobile and immobile zones simultaneously, which are related by first-order mass exchange. Numerical modeling of reactive transport in this kind of multicontinua media is complex and demanding because of the high dimensionality of the problem. In this paper, we establish the governing equations of reactive transport in multicontinuum media incorporating chemical kinetics into the governing equations. We propose a general numerical solution of reactive transport with MRMT by applying direct substitution approach (DSA) based on Newton-Raphson method. The efficiency of the proposed algorithm benefits of the block structure of the system, which allows us to eliminate immobile zones equations and leads to significant savings in CPU time. We test the validity of the developed solution by comparison with other numerical and analytical solutions.

Internal conversion of singlet and triplet states employing numerical DFT/MRCI derivative couplings: Implementation, tests, and application to xanthone.

We present an efficient implementation of nonadiabatic coupling matrix elements (NACMEs) for density functional theory/multireference configuration interaction (DFT/MRCI) wave functions of singlet and triplet multiplicity and an extension of the Vibes program that allows us to determine rate constants for internal conversion (IC) in addition to intersystem crossing (ISC) nonradiative transitions. Following the suggestion of Plasser et al. [J. Chem. Theory Comput. 12, 1207 (2016)], the derivative couplings are computed as finite differences of wave function overlaps. Several measures have been taken to speed up the calculation of the NACMEs. Schur's determinant complement is employed to build up the determinant of the full matrix of spin-blocked orbital overlaps from precomputed spin factors with fixed orbital occupation. Test calculations on formaldehyde, pyrazine, and xanthone show that the mutual excitation level of the configurations at the reference and displaced geometries can be restricted to 1. In combination with a cutoff parameter of tnorm = 10-8 for the DFT/MRCI wave function expansion, this approximation leads to substantial savings of cpu time without essential loss of precision. With regard to applications, the photoexcitation decay kinetics of xanthone in apolar media and in aqueous solution is in the focus of the present work. The results of our computational study substantiate the conjecture that S1 T2 reverse ISC outcompetes the T2 ↝ T1 IC in aqueous solution, thus explaining the occurrence of delayed fluorescence in addition to prompt fluorescence.

A Quasi-Monte Carlo method to compute scattering effects in radiative heat transfer: Application to a sooted jet flame

The Monte-Carlo simulation of radiative heat transfer is known to be very accurate. Its convergence can be significantly improved using Randomized Quasi-Monte-Carlo (RQMC) methods that rely on low-discrepancy sequences. The RQMC approach derived recently to deal with thermal radiation is, however, limited to non-scattering media. The present work proposes a methodology aiming at extending this technique to scattering media based on the prior estimation of the low-discrepancy sequence dimension. Firstly, the method is tested on 3D homogeneous fields with various operating points based on the domain’s optical thickness and albedo. It is observed that, for a given number of generated rays, the error can be reduced by up to one order of magnitude. Secondly, the RQMC approach is combined with importance sampling to increase its efficiency further. The number of rays required is even lower, resulting in saving CPU time to reach a given error. The RQMC approach is then applied along with an accurate model for soot particles’ radiative properties: the Rayleigh-Debye-Gans for Fractal Aggregates (RDGFA) theory. The model assumes a complex morphological shape of particles contrary to Rayleigh theory that is valid for spherical particles only. Monte-Carlo simulations are performed on a fixed turbulent sooted flame field taken from coupled calculations with large-eddy simulation. The overall CPU cost is divided by a factor of 2 compared to a standard Monte Carlo calculation. The simulation allows for accurate quantification of soot scattering effects with RDGFA, which eventually appear small in this configuration. On the other hand, the use of Rayleigh theory strongly underpredicts the actual scattering impact.

Analysis of electromagnetic non-destructive evaluation modelling using Stratton-Chu formulation-based fast algorithm.

The eddy current non-destructive evaluation (NDE) modelling using Stratton-Chu formulation-based fast algorithm is analysed. Stratton-Chu formulations, which have no low frequency breakdown issue, are selected for modelling electromagnetic NDE problems with low frequency and high conductivity approximations. As the main contribution of this article, the robustness and efficiency of the approximations, which result in big savings in both memory and CPU time, are validated and analysed using examples from practical EC testing. The boundary element method (BEM) is used to discretize the integral equations into a linear system of equations: the first order Rao-Wilton-Glisson (RWG) vector basis functions with the flat triangle meshes of the object and pulse basis functions are selected to expand the equivalent surface currents and the normal component of magnetic fields, respectively. Then the multilevel adaptive cross approximation (MLACA) algorithm is applied to accelerate the iterative solution process. The performance and efficiency of adaptively applying a multi-stage (level) algorithm based on the criteria concluded for the operators are shown. This article is part of the theme issue 'Advanced electromagnetic non-destructive evaluation and smart monitoring'.

Analytical Expressions for Lightning Electromagnetic Fields With Arbitrary Channel-Base Current. Part II: Validation and Computational Performance

This article provides the implementation procedure, the validation, and some considerations on the computational efforts of the analytical expressions developed for the lightning electromagnetic fields presented in the companion paper. The validation is presented with different configurations in terms of channel-base current, ground conductivity, and distance to the lightning channel, comparing the obtained results with the numerical integration of the classical formulas. The comparative analysis shows a perfect agreement between the proposed analytical approach and reference numerical simulations. Moreover, the computational effort of the proposed method is discussed, focusing the attention on the choice of the points in which the channel has to be divided in order to maximize the CPU time savings without losing accuracy.

Updating strategy of a domain decomposition preconditioner for parallel solution of dynamic fracture problems

Iterative methods used for analysis of multiphysics problems that result in localized solution features may be highly demanding from a computational standpoint and often require special treatment to be more efficient. Dynamic fracture of metals is one such example in which a nonlinear thermo-mechanical system leads to strain localization, resulting in shear bands and/or cracks, in which iterative solvers may have difficult time converging.In the current paper, we develop a novel updating domain decomposition preconditioner for parallel solution of dynamic fracture problems, in which shearbands are modeled using a temperature-dependent viscoplastic material model and fracture is modeled by the phase-field method. The domain decomposition method is based on the Additive Schwarz Method (ASM).The key idea is to decompose the computational domain into two subdomains, a localized subdomain that includes all localized features of the solution and a healthy subdomain for the remaining part of the domain. In this way, one can apply different solvers in each subdomain, i.e. focus more effort in the localized subdomain. In this work, an LU solver is applied in both subdomains, however, while the localized subdomain is solved exactly at every nonlinear iteration, the healthy subdomain LU operator is reused and only selectively updated. Hence, significant CPU time savings associated with the setup of the preconditioner can be achieved.To this end, we propose a strategy for updating the preconditioner in the healthy subdomain. The strategy is based on an idealized performance-based optimization procedure that takes into account machine on-the-fly execution time. Three dynamic fracture problems corresponding to different failure modes are investigated. Excellent performance of the proposed updating preconditioner is reported in serial and parallel computers.

Design of Ku- and Ka-Band Flat Dual Circular Polarized Reflectarrays by Combining Variable Rotation Technique and Element Size Variation

In this work circularly polarized flat reflectarray antennas are designed to work at downlink frequencies in Ku- and Ka-band. The design approach for circular polarization radiation combines both the variable rotation technique (VRT) and the element size variation technique (ESVT). The reflectarray element employed consists of a rotable split ring operating at 12 GHz (Ku-band) that surrounds two orthogonal sets of parallel dipoles of variable size operating at 20 GHz (Ka-band). Rotations of split rings and different sizes of dipoles are exploited as degrees of freedom by VRT and ESVT for designing focused beams dual band dual circular polarized reflectarrays. The split ring can be used to achieve single focused beams for circular polarization by VRT at 12 GHz, and the set of orthogonal dipoles can be used to achieve two focused beams for dual circular polarization by ESVT at 20 GHz. In this way, the reflectarray element has been used to design two flat dual band circularly polarized reflectarrays that generate three pencil beams either in the same direction or in different directions: one pencil beam with right-hand circular polarization (RHCP) at 12 GHz and two pencil beams with both RHCP and left-hand circular polarization (LHCP) at 20 GHz. The designed reflectarrays are carried out under local periodicity assumption by means of in-house electromagnetic software that applies the Method of Moments in the Spectral Domain. Validations of the in-house electromagnetic software show significant CPU time savings with respect to CPU time consumption provided by commercial software of general purpose as CST electromagnetics software. Numerical results of radiation patterns produced by the designed reflectarrays show efficiencies around 65% and a bandwidth of 6% for a main beam cross-polarization discrimination of 25 dB. The numerical results are satisfactorily validated by CST software and a tolerance errors study has been numerically carried out with acceptable results.

A high-resolution method based on off-step non-polynomial spline approximations for the solution of Burgers-Fisher and coupled nonlinear Burgers’ equations

PurposeThis paper aims to develop a new high accuracy numerical method based on off-step non-polynomial spline in tension approximations for the solution of Burgers-Fisher and coupled nonlinear Burgers’ equations on a graded mesh. The spline method reported here is third order accurate in space and second order accurate in time. The proposed spline method involves only two off-step points and a central point on a graded mesh. The method is two-level implicit in nature and directly derived from the continuity condition of the first order space derivative of the non-polynomial tension spline function. The linear stability analysis of the proposed method has been examined and it is shown that the proposed two-level method is unconditionally stable for a linear model problem. The method is directly applicable to problems in polar systems. To demonstrate the strength and utility of the proposed method, the authors have solved the generalized Burgers-Huxley equation, generalized Burgers-Fisher equation, coupled Burgers-equations and parabolic equation in polar coordinates. The authors show that the proposed method enables us to obtain the high accurate solution for high Reynolds number.Design/methodology/approachIn this method, the authors use only two-level in time-direction, and at each time-level, the authors use three grid points for the unknown function u(x,t) and two off-step points for the known variable x in spatial direction. The methodology followed in this paper is the construction of a non-polynomial spline function and using its continuity properties to obtain consistency condition, which is third order accurate on a graded mesh and fourth order accurate on a uniform mesh. From this consistency condition, the authors derive the proposed numerical method. The proposed method, when applied to a linear equation is shown to be unconditionally stable. To assess the validity and accuracy, the method is applied to solve several benchmark problems, and numerical results are provided to demonstrate the usefulness of the proposed method.FindingsThe paper provides a third order numerical scheme on a graded mesh and fourth order spline method on a uniform mesh obtained directly from the consistency condition. In earlier methods, consistency conditions were only second order accurate. This brings an edge over other past methods. Also, the method is directly applicable to physical problems involving singular coefficients. So no modification in the method is required at singular points. This saves CPU time and computational costs.Research limitations/implicationsThere are no limitations. Obtaining a high accuracy spline method directly from the consistency condition is a new work. Also being an implicit method, this method is unconditionally stable.Practical implicationsPhysical problems with singular and non-singular coefficients are directly solved by this method.Originality/valueThe paper develops a new method based on non-polynomial spline approximations of order two in time and three (four) in space, which is original and has lot of value because many benchmark problems of physical significance are solved in this method.