Note on derivation of order conditions for ARKN methods for perturbed oscillators

J.M. Franco derived the sufficient order conditions as well as the necessary and sufficient order conditions for his Adapted Runge-Kutta-Nyström methods (in short notation ARKN methods) based on the B-series theory [J.M. Franco, Runge-Kutta-Nyström methods adapted to the numerical integration of perturbed oscillators, Comput. Phys. Comm. 147 (2002) 770-787]. Unfortunately, some fundamental mistakes have been made in the derivation of order conditions in that paper. In view of importance of the algebraic order theory for ARKN methods, a new and correct derivation of the order conditions for the ARKN methods is presented in this short note.     

ERKN integrators for systems of oscillatory second-order differential equations

For systems of oscillatory second-order differential equations y^″+My=f with M∈Rm×m, a symmetric positive semi-definite matrix, X. Wu et al. have proposed the multidimensional ARKN methods [X. Wu, X. You, J. Xia, Order conditions for ARKN methods solving oscillatory systems, Comput. Phys. Comm. 180 (2009) 2250-2257], which are an essential generalization of J.M. Franco's ARKN methods for one-dimensional problems or for systems with a diagonal matrix M=w²I [J.M. Franco, Runge-Kutta-Nyström methods adapted to the numerical integration of perturbed oscillators, Comput. Phys. Comm. 147 (2002) 770-787]. One of the merits of these methods is that they integrate exactly the unperturbed oscillators y^″+My=0. Regretfully, even for the unperturbed oscillators the internal stages Yi of an ARKN method fail to equal the values of the exact solution y(t) at tn+cih, respectively. Recently H. Yang et al. proposed the ERKN methods to overcome this drawback [H.L. Yang, X.Y. Wu, Xiong You, Yonglei Fang, Extended RKN-type methods for numerical integration of perturbed oscillators, Comput. Phys. Comm. 180 (2009) 1777-1794]. However, the ERKN methods in that paper are only considered for the special case where M is a diagonal matrix with nonnegative entries. The purpose of this paper is to extend the ERKN methods to the general case with M∈Rm×m, and the perturbing function f depends only on y. Numerical experiments accompanied demonstrates that the ERKN methods are more efficient than the existing methods for the computation of oscillatory systems. In particular, if M∈Rm×m is a symmetric positive semi-definite matrix, it is highly important for the new ERKN integrators to show the energy conservation in the numerical experiments for problems with Hamiltonian in comparison with the well-known methods in the scientific literature. Those so called separable Hamiltonians arise in many areas of physical sciences, e.g., macromolecular dynamics, astronomy, and classical mechanics.     

Order conditions for ARKN methods solving oscillatory systems

For the perturbed oscillators in one-dimensional case, J.M. Franco designed the so-called Adapted Runge-Kutta-Nyström (ARKN) methods and derived the sufficient order conditions as well as the necessary and sufficient order conditions for ARKN methods based on the B-series theory [J.M. Franco, Runge-Kutta-Nyström methods adapted to the numerical integration of perturbed oscillators, Comput. Phys. Comm. 147 (2002) 770-787]. These methods integrate exactly the unperturbed oscillators and are highly efficient when the perturbing function is small. Unfortunately, some critical mistakes have been made in the derivation of order conditions in that paper. On the basis of the results from that paper, Franco extended directly the ARKN methods and the corresponding order conditions to multidimensional case where the perturbed function f does not depend on the first derivative y^′ [J.M. Franco, New methods for oscillatory systems based on ARKN methods, Appl. Numer. Math. 56 (2006) 1040-1053]. In this paper, we present the order conditions for the ARKN methods for the general multidimensional perturbed oscillators where the perturbed function f may depend on only y or on both y and y^′.     

A DAFT DL_POLY distributed memory adaptation of the Smoothed Particle Mesh Ewald method

The Smoothed Particle Mesh Ewald method [U. Essmann, L. Perera, M.L. Berkowtz, T. Darden, H. Lee, L.G. Pedersen, J. Chem. Phys. 103 (1995) 8577] for calculating long ranged forces in molecular simulation has been adapted for the parallel molecular dynamics code DL_POLY_3 [I.T. Todorov, W. Smith, Philos. Trans. Roy. Soc. London 362 (2004) 1835], making use of a novel 3D Fast Fourier Transform (DAFT) [I.J. Bush, The Daresbury Advanced Fourier transform, Daresbury Laboratory, 1999] that perfectly matches the Domain Decomposition (DD) parallelisation strategy [W. Smith, Comput. Phys. Comm. 62 (1991) 229; M.R.S. Pinches, D. Tildesley, W. Smith, Mol. Sim. 6 (1991) 51; D. Rapaport, Comput. Phys. Comm. 62 (1991) 217] of the DL_POLY_3 code. In this article we describe software adaptations undertaken to import this functionality and provide a review of its performance.     

Coupled cluster algorithms for networks of shared memory parallel processors

As the popularity of using SMP systems as the building blocks for high performance supercomputers increases, so too increases the need for applications that can utilize the multiple levels of parallelism available in clusters of SMPs. This paper presents a dual-layer distributed algorithm, using both shared-memory and distributed-memory techniques to parallelize a very important algorithm (often called the “gold standard”) used in computational chemistry, the single and double excitation coupled cluster method with perturbative triples, i.e. CCSD(T). The algorithm is presented within the framework of the GAMESS [M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.J. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen, S. Su, T.L. Windus, M. Dupuis, J.A. Montgomery, General atomic and molecular electronic structure system, J. Comput. Chem. 14 (1993) 1347–1363]. (General Atomic and Molecular Electronic Structure System) program suite and the Distributed Data Interface [M.W. Schmidt, G.D. Fletcher, B.M. Bode, M.S. Gordon, The distributed data interface in GAMESS, Comput. Phys. Comm. 128 (2000) 190]. (DDI), however, the essential features of the algorithm (data distribution, load-balancing and communication overhead) can be applied to more general computational problems. Timing and performance data for our dual-level algorithm is presented on several large-scale clusters of SMPs.     

Spatially adaptive techniques for level set methods and incompressible flow

Since the seminal work of [Sussman, M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59] on coupling the level set method of [Osher S, Sethian J. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79:12–49] to the equations for two-phase incompressible flow, there has been a great deal of interest in this area. That work demonstrated the most powerful aspects of the level set method, i.e. automatic handling of topological changes such as merging and pinching, as well as robust geometric information such as normals and curvature. Interestingly, this work also demonstrated the largest weakness of the level set method, i.e. mass or information loss characteristic of most Eulerian capturing techniques. In fact, [Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59] introduced a partial differential equation for battling this weakness, without which their work would not have been possible. In this paper, we discuss both historical and most recent works focused on improving the computational accuracy of the level set method focusing in part on applications related to incompressible flow due to both of its popularity and stringent accuracy requirements. Thus, we discuss higher order accurate numerical methods such as Hamilton–Jacobi WENO [Jiang G-S, Peng D. Weighted ENO schemes for Hamilton–Jacobi equations. SIAM J Sci Comput 2000;21:2126–43], methods for maintaining a signed distance function, hybrid methods such as the particle level set method [Enright D, Fedkiw R, Ferziger J, Mitchell I. A hybrid particle level set method for improved interface capturing. J Comput Phys 2002;183:83–116] and the coupled level set volume of fluid method [Sussman M, Puckett EG. A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows. J Comput Phys 2000;162:301–37], and adaptive gridding techniques such as the octree approach to free surface flows proposed in [Losasso F, Gibou F, Fedkiw R. Simulating water and smoke with an octree data structure, ACM Trans Graph (SIGGRAPH Proc) 2004;23:457–62].     

An analysis of and a comparison between the discontinuous Galerkin and the spectral finite volume methods

The discontinuous Galerkin method has been developed and applied extensively to solve hyperbolic conservation laws in recent years. More recently Wang et al. developed a class of discontinuous Petrov–Galerkin method, termed spectral (finite) volume method [J. Comput. Phys. 78 (2002) 210; J. Comput. Phys. 179 (2002) 665; J. Sci. Comput. 20 (2004) 137]. In this paper we perform a Fourier type analysis on both methods when solving linear one-dimensional conservation laws. A comparison between the two methods is given in terms of accuracy, stability, and convergence. Numerical experiments are performed to validate this analysis and comparison.     

Extended RKN methods with FSAL property for oscillatory systems

In this paper, we present extended RKN methods with FSAL property for numerical integration of perturbed oscillators. These numerical integrators are based on the extended Runge-Kutta-Nyström-type methods proposed by Yang et al. [H.L. Yang, X.Y. Wu, Xiong You, Y.L. Fang, Extended RKN-type methods for numerical integration of perturbed oscillators, Comput. Phys. Comm. 180 (2009) 1777-1794]. The numerical stability and phase property of the new methods are analyzed. The paper is accompanied by numerical experiments that demonstrate the effectiveness and robustness of our new methods in comparison with some well-known methods appeared in the scientific literature.     

Addendum to: “The cross sections for one phonon emission and absorption by slow neutrons in superfluid liquid Helium” [Comput. Phys. Comm. 151 (2) (2003) 141-148]: Angular distribution of UCN

The UCN differential angular distribution produced by superfluid liquid helium down-scattering was numerically computed in conjunction with the previous works of the total cross sections [H. Yoshiki, Comput. Phys. Comm. 151 (2003) 141-148]. The distribution for the maximum UCN energy limit of 252 neV (Be barrier potential) is very nearly isotropic.     

ExHuME 1.3: A Monte Carlo event generator for exclusive diffraction

We have written the Exclusive Hadronic Monte Carlo Event (ExHuME) generator. ExHuME is based around the perturbative QCD calculation of Khoze, Martin and Ryskin of the process pp→p+X+p, where X is a centrally produced colour singlet system.

Program summary

Title of program:ExHuMECatalogue identifier:ADYA_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADYA_v1_0Program obtainable from:CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions:NoneProgramming language used:C++, some FORTRANComputer:Any computer with UNIX capability. Users should refer to the README file distributed with the source code for further detailsOperating system:Linux, Mac OS XNo. of lines in distributed program, including test data, etc.:111 145No. of bytes in distributed program, including test data, etc.: 791 085Distribution format:tar.gzRAM:60 MBExternal routines/libraries:LHAPDF [http://durpdg.dur.ac.uk/lhapdf/], CLHEP v1.8 or v1.9 [L. Lönnblad, Comput. Phys. Comm. 84 (1994) 307; http://wwwinfo.cern.ch/asd/lhc++/clhep/]Subprograms used:Pythia [T. Sjostrand et al., Comput. Phys. Comm. 135 (2001) 238], HDECAY [A. Djouadi, J. Kalinowski, M. Spira, HDECAY: A program for Higgs boson decays in the standard model and its supersymmetric extension, Comput. Phys. Comm. 108 (1998) 56, hep-ph/9704448]. Both are distributed with the source codeNature of problem:Central exclusive production offers the opportunity to study particle production in a uniquely clean environment for a hadron collider. This program implements the KMR model [V.A. Khoze, A.D. Martin, M.G. Ryskin, Prospects for New Physics observations in diffractive processes at the LHC and Tevatron, Eur. Phys. J. C 23 (2002) 311, hep-ph/0111078], which is the only fully perturbative model of exclusive production.Solution method:Monte Carlo techniques are used to produce the central exclusive parton level system. Pythia routines are then used to develop a realistic hadronic system.Restrictions:The program is, at present, limited to Higgs, di-gluon and di-quark production. However, in principle it is not difficult to include more.Running time:Approximately 10 minutes for 10000 Higgs events on an Apple 1 GHz G4 PowerPC.     

Computation of incompressible bubble dynamics with a stabilized finite element level set method

This paper presents a stabilized finite element method for the three dimensional computation of incompressible bubble dynamics using a level set method. The interface between the two phases is resolved using the level set approach developed by Sethian [Level Set Methods and Fast Marching Methods, Cambridge University Press, 1999], Sussman et al. [J. Comput. Phys. 114 (1994) 146], and Sussman et al. [J. Comput. Phys. 148 (1999) 81–124]. In this approach the interface is represented as a zero level set of a smooth function. The streamline-upwind/Petrov–Galerkin method was used to discretize the governing flow and level set equations. The continuum surface force (CSF) model proposed by Brackbill et al. [J. Comput. Phys. 100 (1992) 335–354] was applied in order to account for surface tension effects. To restrict the interface from moving while re-distancing, an improved re-distancing scheme proposed in the finite difference context [J. Comput. Phys. 148 (1999) 81–124] is adapted for finite element discretization. This enables us to accurately compute the flows with large density and viscosity differences, as well as surface tension. The capability of the resultant algorithm is demonstrated with two and three dimensional numerical examples of a single bubble rising through a quiescent liquid, and two bubble coalescence.     

Numerical solution of the shallow water equations with a fractional step method

A fractional step technique for the numerical solution of the shallow water equations is applied to study the evolution of the potential vorticity field. The height and velocity field of the shallow water equations are discretized on a fixed Eulerian grid and time-stepped with a fractional step method recently reported in [M. Shoucri, Comput. Phys. Comm. 164 (2004) 396; M. Shoucri, A. Qaddouri, M. Tanguay, J. Côté, A Fractional Steps Method for the Numerical Solution of the Shallow Water Equations, International Workshop on Solution of Partial Differential Equations, The Fields Institute, Toronto, August 2002], where the Riemann invariants of the equations are interpolated at each time step along the characteristics using a cubic spline interpolation. The potential vorticity, which develops steep gradients and evolve into thin filaments during the evolution, is nicely calculated at every time-step from the solution of the code. The method is efficient and has lower numerical diffusion than other methods, since it evolves the equations without the iterative steps involved in the multi-dimensional interpolation problem, and without the iteration associated with the intermediate step of solving a Helmholtz equation, usually associated with other methods like the semi-Lagrangian method. The absence of iterative steps in the present technique makes it very suitable for problems in which small time steps and grid sizes are required, as for instance in the present problem where steepness of the gradients and small scale structures are the main features of the potential vorticity, and more generally for problems of regional climate modeling. The simplicity of the method makes it very suitable for parallel computer.     

A comparable study of image approximations to the reaction field

The recently developed high-order accurate multiple image approximation to the reaction field for a charge inside a dielectric sphere [J. Comput. Phys. 223 (2007) 846-864] is compared favorably to other commonly employed reaction field schemes. These methods are of particular interest because they are useful in the study of biological macromolecules by the Monte Carlo and Molecular Dynamics methods.     

Extended RKN-type methods with minimal dispersion error for perturbed oscillators

In this paper, extended Runge-Kutta-Nyström-type methods with minimal dispersion error for the numerical integration of perturbed oscillators are presented, which are based on the order conditions derived by Yang et al. (H.L. Yang, X.Y. Wu, X. You, Y.L. Fang, Extended RKN-type methods for numerical integration of perturbed oscillators, Comput. Phys. Comm. 180 (2009) 1777-1794). The numerical stability and phase properties of the new methods are analyzed. Numerical experiments are reported in comparison with some well-known high quality codes proposed in the scientific literature.     

Sixth-order symmetric and symplectic exponentially fitted modified Runge-Kutta methods of Gauss type

The construction of symmetric and symplectic exponentially fitted modified Runge-Kutta (RK) methods for the numerical integration of Hamiltonian systems with oscillatory solutions is considered. In a previous paper [H. Van de Vyver, A fourth order symplectic exponentially fitted integrator, Comput. Phys. Comm. 176 (2006) 255-262] a two-stage fourth-order symplectic exponentially fitted modified RK method has been proposed. Here, two three-stage symmetric and symplectic exponentially fitted integrators of Gauss type, either with fixed nodes or variable nodes, are derived. The algebraic order of the new integrators is also analyzed, obtaining that they possess sixth-order as the classical three-stage RK Gauss method. Numerical experiments with some oscillatory problems are presented to show that the new methods are more efficient than other symplectic RK Gauss codes proposed in the scientific literature.     

DPEMC: A Monte Carlo for double diffraction

We extend the POMWIG Monte Carlo generator developed by B. Cox and J. Forshaw, to include new models of central production through inclusive and exclusive double Pomeron exchange in proton-proton collisions. Double photon exchange processes are described as well, both in proton-proton and heavy-ion collisions. In all contexts, various models have been implemented, allowing for comparisons and uncertainty evaluation and enabling detailed experimental simulations.

Program summary

Title of the program:DPEMC, version 2.4Catalogue identifier: ADVFProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADVFProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandComputer: any computer with the FORTRAN 77 compiler under the UNIX or Linux operating systemsOperating system: UNIX; LinuxProgramming language used: FORTRAN 77High speed storage required:<25 MBNo. of lines in distributed program, including test data, etc.: 71 399No. of bytes in distributed program, including test data, etc.: 639 950Distribution format: tar.gzNature of the physical problem: Proton diffraction at hadron colliders can manifest itself in many forms, and a variety of models exist that attempt to describe it [A. Bialas, P.V. Landshoff, Phys. Lett. B 256 (1991) 540; A. Bialas, W. Szeremeta, Phys. Lett. B 296 (1992) 191; A. Bialas, R.A. Janik, Z. Phys. C 62 (1994) 487; M. Boonekamp, R. Peschanski, C. Royon, Phys. Rev. Lett. 87 (2001) 251806; Nucl. Phys. B 669 (2003) 277; R. Enberg, G. Ingelman, A. Kissavos, N. Timneanu, Phys. Rev. Lett. 89 (2002) 081801; R. Enberg, G. Ingelman, L. Motyka, Phys. Lett. B 524 (2002) 273; R. Enberg, G. Ingelman, N. Timneanu, Phys. Rev. D 67 (2003) 011301; B. Cox, J. Forshaw, Comput. Phys. Comm. 144 (2002) 104; B. Cox, J. Forshaw, B. Heinemann, Phys. Lett. B 540 (2002) 26; V. Khoze, A. Martin, M. Ryskin, Phys. Lett. B 401 (1997) 330; Eur. Phys. J. C 14 (2000) 525; Eur. Phys. J. C 19 (2001) 477; Erratum, Eur. Phys. J. C 20 (2001) 599; Eur. Phys. J. C 23 (2002) 311]. This program implements some of the more significant ones, enabling the simulation of central particle production through color singlet exchange between interacting protons or antiprotons.Method of solution: The Monte Carlo method is used to simulate all elementary 2→2 and 2→1 processes available in HERWIG. The color singlet exchanges implemented in DPEMC are implemented as functions reweighting the photon flux already present in HERWIG.Restriction on the complexity of the problem: The program relying extensively on HERWIG, the limitations are the same as in [G. Marchesini, B.R. Webber, G. Abbiendi, I.G. Knowles, M.H. Seymour, L. Stanco, Comput. Phys. Comm. 67 (1992) 465; G. Corcella, I.G. Knowles, G. Marchesini, S. Moretti, K. Odagiri, P. Richardson, M. Seymour, B. Webber, JHEP 0101 (2001) 010].Typical running time: Approximate times on a 800 MHz Pentium III: 5-20 min per 10 000 unweighted events, depending on the process under consideration.     

On improving mass-conservation properties of the hybrid particle-level-set method

For the computation of multi-phase flows level-set methods are an attractive alternative to volume-of-fluid or front-tracking approaches. For improving their accuracy and efficiency the hybrid particle-level-set modification was proposed by Enright et al. [Enright D, Fedkiw R, Ferziger J, Mitchell I. A hybrid particle-level-set method for improved interface capturing. J Comput Phys 2002;183:83-116]. In actual applications the overall properties of a level-set method, such as mass conservation, are strongly affected by discretization schemes and algorithmic details. In this paper we address these issues with the objective of determining the optimum alternatives for the purpose of direct numerical simulation of dispersed-droplet flows. We evaluate different discretization schemes for curvature and unit normal vector at the interface. Another issue is the particular formulation of the reinitialization of the level-set function which significantly affects the quality of computational results. Different approaches employing higher-order schemes for discretization, supplemented either by a correction step using marker particles (Enright et al., 2002) or by additional constraints [Sussman M, Almgren AS, Bell JB, Colella P, Howell LH, Welcome ML. An adaptive level set approach for incompressible two-phase flows. J Comput Phys 1999;148:81-124] are analyzed. Different parameter choices for the hybrid particle-level-set method are evaluated with the purpose of increasing the efficiency of the method. Aiming at large-scale computations we find that in comparison with pure level-set methods the hybrid particle-level-set method exhibits better mass-conservation properties, especially in the case of marginally resolved interfaces.     

ADER schemes on unstructured meshes for nonconservative hyperbolic systems: Applications to geophysical flows

We develop a new family of well-balanced path-conservative quadrature-free one-step ADER finite volume and discontinuous Galerkin finite element schemes on unstructured meshes for the solution of hyperbolic partial differential equations with non-conservative products and stiff source terms. The fully discrete formulation is derived using the recently developed framework of explicit one-step P_NP_M schemes of arbitrary high order of accuracy in space and time for conservative hyperbolic systems [Dumbser M, Balsara D, Toro EF, Munz CD. A unified framework for the construction of one-step finite-volume and discontinuous Galerkin schemes. J Comput Phys 2008;227:8209–53]. The two key ingredients of our high order approach are: first, the high order accurate P_NP_M reconstruction operator on unstructured meshes, using the WENO strategy presented in [Dumbser M, Käser M, Titarev VA Toro EF. Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinear hyperbolic systems. J Comput Phys 2007;226:204–43] to ensure monotonicity at discontinuities, and second, a local space–time Galerkin scheme to predict the evolution of the reconstructed polynomial data inside each element during one time step to obtain a high order accurate one-step time discretization. This approach is also able to deal with stiff source terms as shown in [Dumbser M, Enaux C, Toro EF. Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws. J Comput Phys 2008;227:3971–4001]. These two key ingredients are combined with the recently developed path-conservative methods of Parés [Parés C. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM J Numer Anal 2006;44:300–21] and Castro et al. [Castro MJ, Gallardo JM, Parés C. High-order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. Applications to shallow-water systems. Math Comput 2006;75:1103–34] to treat the non-conservative products properly. We show applications of our method to the two-layer shallow water equations as well as applications to the recently published depth-averaged two-fluid flow model of Pitman and Le [Pitman EB, Le L. A two-fluid model for avalanche and debris flows. Philos Trans Roy Soc A 2005;363:1573–601].     

Collective mode formulation of the response algorithm for solving Kohn-Sham equations

We have developed a “collective mode” version of the response-iteration algorithm for solving the nonlinear Kohn-Sham equations. The algorithm utilizes approximation methods for the density-density response function that are known from microscopic many-body theories of strongly interacting Fermi systems. The major advantage over our previously proposed algorithm (J. Auer and E. Krotscheck, Comput. Phys. Comm. 118 (1999) 139-144) is that the new method needs the computation of occupied states only. Using spherical jellium clusters with up to 2000 electrons as an example, we show that the approximations implicit to our new algorithms do not deteriorate the convergence rate. An even simpler version approximates the density-density response function by that of a charged bose gas with the same density. This algorithm converges somewhat more slowly, but still provides a viable method for solving Kohn-Sham equations for small clusters.