@article{Reich2006, author = {Reich, Sebastian}, title = {Linearly implicit time stepping methods for numerical weather prediction}, series = {BIT : numerical mathematics ; the leading applied mathematics journal for all computational mathematicians}, volume = {46}, journal = {BIT : numerical mathematics ; the leading applied mathematics journal for all computational mathematicians}, publisher = {Springer}, address = {Dordrecht}, issn = {0006-3835}, doi = {10.1007/s10543-006-0065-0}, pages = {607 -- 616}, year = {2006}, abstract = {The efficient time integration of the dynamic core equations for numerical weather prediction (NWP) remains a key challenge. One of the most popular methods is currently provided by implementations of the semi-implicit semi-Lagrangian (SISL) method, originally proposed by Robert (J. Meteorol. Soc. Jpn., 1982). Practical implementations of the SISL method are, however, not without certain shortcomings with regard to accuracy, conservation properties and stability. Based on recent work by Gottwald, Frank and Reich (LNCSE, Springer, 2002), Frank, Reich, Staniforth, White and Wood (Atm. Sci. Lett., 2005) and Wood, Staniforth and Reich (Atm. Sci. Lett., 2006) we propose an alternative semi-Lagrangian implementation based on a set of regularized equations and the popular Stormer-Verlet time stepping method in the context of the shallow-water equations (SWEs). Ultimately, the goal is to develop practical implementations for the 3D Euler equations that overcome some or all shortcomings of current SISL implementations.}, language = {en} } @article{FrankMooreReich2006, author = {Frank, Jason and Moore, Brian E. and Reich, Sebastian}, title = {Linear PDEs and numerical methods that preserve a multisymplectic conservation law}, issn = {1064-8275}, doi = {10.1137/050628271}, year = {2006}, abstract = {Multisymplectic methods have recently been proposed as a generalization of symplectic ODE methods to the case of Hamiltonian PDEs. Their excellent long time behavior for a variety of Hamiltonian wave equations has been demonstrated in a number of numerical studies. A theoretical investigation and justification of multisymplectic methods is still largely missing. In this paper, we study linear multisymplectic PDEs and their discretization by means of numerical dispersion relations. It is found that multisymplectic methods in the sense of Bridges and Reich [Phys. Lett. A, 284 ( 2001), pp. 184-193] and Reich [J. Comput. Phys., 157 (2000), pp. 473-499], such as Gauss-Legendre Runge-Kutta methods, possess a number of desirable properties such as nonexistence of spurious roots and conservation of the sign of the group velocity. A certain CFL-type restriction on Delta t/Delta x might be required for methods higher than second order in time. It is also demonstrated by means of the explicit midpoint method that multistep methods may exhibit spurious roots in the numerical dispersion relation for any value of Delta t/Delta x despite being multisymplectic in the sense of discrete variational mechanics [J. E. Marsden, G. P. Patrick, and S. Shkoller, Commun. Math. Phys., 199 (1999), pp. 351-395]}, language = {en} } @article{BridgesReich2006, author = {Bridges, Thomas J. and Reich, Sebastian}, title = {Numerical methods for Hamiltonian PDEs}, issn = {0305-4470}, doi = {10.1088/0305-4470/39/19/S02}, year = {2006}, abstract = {The paper provides an introduction and survey of conservative discretization methods for Hamiltonian partial differential equations. The emphasis is on variational, symplectic and multi-symplectic methods. The derivation of methods as well as some of their fundamental geometric properties are discussed. Basic principles are illustrated by means of examples from wave and fluid dynamics}, language = {en} } @article{AkhmatskayaBouRabeeReich2009, author = {Akhmatskaya, Elena and Bou-Rabee, Nawaf and Reich, Sebastian}, title = {A comparison of generalized hybrid Monte Carlo methods with and without momentum flip}, issn = {0021-9991}, doi = {10.1016/j.jcp.2008.12.014}, year = {2009}, abstract = {The generalized hybrid Monte Carlo (GHMC) method combines Metropolis corrected constant energy simulations with a partial random refreshment step in the particle momenta. The standard detailed balance condition requires that momenta are negated upon rejection of a molecular dynamics proposal step. The implication is a trajectory reversal upon rejection, which is undesirable when interpreting GHMC as thermostated molecular dynamics. We show that a modified detailed balance condition can be used to implement GHMC without momentum flips. The same modification can be applied to the generalized shadow hybrid Monte Carlo (GSHMC) method. Numerical results indicate that GHMC/GSHMC implementations with momentum flip display a favorable behavior in terms of sampling efficiency, i.e., the traditional GHMC/GSHMC implementations with momentum flip got the advantage of a higher acceptance rate and faster decorrelation of Monte Carlo samples. The difference is more pronounced for GHMC. We also numerically investigate the behavior of the GHMC method as a Langevin-type thermostat. We find that the GHMC method without momentum flip interferes less with the underlying stochastic molecular dynamics in terms of autocorrelation functions and it to be preferred over the GHMC method with momentum flip. The same finding applies to GSHMC.}, language = {en} } @article{AkhmatskayaBouRabeeReich2009, author = {Akhmatskaya, Elena and Bou-Rabee, Nawaf and Reich, Sebastian}, title = {Erratum to "A comparison of generalized hybrid Monte Carlo methods with and without momentum flip" [J. Comput. Phys. 228 (2009), S. 2256 - 2265]}, issn = {0021-9991}, doi = {10.1016/j.jcp.2009.06.039}, year = {2009}, abstract = {The generalized hybrid Monte Carlo (GHMC) method combines Metropolis corrected constant energy simulations with a partial random refreshment step in the particle momenta. The standard detailed balance condition requires that momenta are negated upon rejection of a molecular dynamics proposal step. The implication is a trajectory reversal upon rejection, which is undesirable when interpreting GHMC as thermostated molecular dynamics. We show that a modified detailed balance condition can be used to implement GHMC without momentum flips. The same modification can be applied to the generalized shadow hybrid Monte Carlo (GSHMC) method. Numerical results indicate that GHMC/GSHMC implementations with momentum flip display a favorable behavior in terms of sampling efficiency, i.e., the traditional GHMC/GSHMC implementations with momentum flip got the advantage of a higher acceptance rate and faster decorrelation of Monte Carlo samples. The difference is more pronounced for GHMC. We also numerically investigate the behavior of the GHMC method as a Langevin-type thermostat. We find that the GHMC method without momentum flip interferes less with the underlying stochastic molecular dynamics in terms of autocorrelation functions and it to be preferred over the GHMC method with momentum flip. The same finding applies to GSHMC.}, language = {en} } @article{LeimkuhlerReich2009, author = {Leimkuhler, Benedict and Reich, Sebastian}, title = {A metropolis adjusted Nos{\´e}-Hoover thermostat}, issn = {0764-583X}, doi = {10.1051/M2an/2009023}, year = {2009}, abstract = {We present a Monte Carlo technique for sampling from the canonical distribution in molecular dynamics. The method is built upon the Nose-Hoover constant temperature formulation and the generalized hybrid Monte Carlo method. In contrast to standard hybrid Monte Carlo methods only the thermostat degree of freedom is stochastically resampled during a Monte Carlo step.}, language = {en} } @article{BergemannGottwaldReich2009, author = {Bergemann, Kay and Gottwald, Georg and Reich, Sebastian}, title = {Ensemble propagation and continuous matrix factorization algorithms}, issn = {0035-9009}, doi = {10.1002/qj.457}, year = {2009}, abstract = {We consider the problem of propagating an ensemble of solutions and its characterization in terms of its mean and covariance matrix. We propose differential equations that lead to a continuous matrix factorization of the ensemble into a generalized singular value decomposition (SVD). The continuous factorization is applied to ensemble propagation under periodic rescaling (ensemble breeding) and under periodic Kalman analysis steps (ensemble Kalman filter). We also use the continuous matrix factorization to perform a re-orthogonalization of the ensemble after each time-step and apply the resulting modified ensemble propagation algorithm to the ensemble Kalman filter. Results from the Lorenz-96 model indicate that the re-orthogonalization of the ensembles leads to improved filter performance.}, language = {en} } @article{CotterHamPainetal.2009, author = {Cotter, Colin J. and Ham, David A. and Pain, Christopher C. and Reich, Sebastian}, title = {LBB stability of a mixed Galerkin finite element pair for fluid flow simulations}, issn = {0021-9991}, doi = {10.1016/j.jcp.2008.09.014}, year = {2009}, abstract = {We introduce a new mixed finite element for solving the 2- and 3-dimensional wave equations and equations of incompressible flow. The element, which we refer to as P1(D)-P2, uses discontinuous piecewise linear functions for velocity and continuous piecewise quadratic functions for pressure. The aim of introducing the mixed formulation is to produce a new flexible element choice for triangular and tetrahedral meshes which satisfies the LBB stability condition and hence has no spurious zero-energy modes. The advantage of this particular element choice is that the mass matrix for velocity is block diagonal so it can be trivially inverted; it also allows the order of the pressure to be increased to quadratic whilst maintaining LBB stability which has benefits in geophysical applications with Coriolis forces. We give a normal mode analysis of the semi-discrete wave equation in one dimension which shows that the element pair is stable, and demonstrate that the element is stable with numerical integrations of the wave equation in two dimensions, an analysis of the resultant discrete Laplace operator in two and three dimensions on various meshes which shows that the element pair does not have any spurious modes. We provide convergence tests for the element pair which confirm that the element is stable since the convergence rate of the numerical solution is quadratic.}, language = {en} } @article{ShinSommerReichetal.2010, author = {Shin, Seoleun and Sommer, Matthias and Reich, Sebastian and N{\´e}vir, Peter}, title = {Evaluation of three spatial discretization schemes with the Galewsky et al. test}, issn = {1530-261X}, doi = {10.1002/Asl.279}, year = {2010}, abstract = {We evaluate the Hamiltonian particle methods (HPM) and the Nambu discretization applied to shallow-water equations on the sphere using the test suggested by Galewsky et al. (2004). Both simulations show excellent conservation of energy and are stable in long-term simulation. We repeat the test also using the ICOSWP scheme to compare with the two conservative spatial discretization schemes. The HPM simulation captures the main features of the reference solution, but wave 5 pattern is dominant in the simulations applied on the ICON grid with relatively low spatial resolutions. Nevertheless, agreement in statistics between the three schemes indicates their qualitatively similar behaviors in the long-term integration.}, language = {en} } @article{BergemannReich2010, author = {Bergemann, Kay and Reich, Sebastian}, title = {A localization technique for ensemble Kalman filters}, issn = {0035-9009}, doi = {10.1002/Qj.591}, year = {2010}, abstract = {Ensemble Kalman filter techniques are widely used to assimilate observations into dynamical models. The phase- space dimension is typically much larger than the number of ensemble members, which leads to inaccurate results in the computed covariance matrices. These inaccuracies can lead, among other things, to spurious long-range correlations, which can be eliminated by Schur-product-based localization techniques. In this article, we propose a new technique for implementing such localization techniques within the class of ensemble transform/square-root Kalman filters. Our approach relies on a continuous embedding of the Kalman filter update for the ensemble members, i.e. we state an ordinary differential equation (ODE) with solutions that, over a unit time interval, are equivalent to the Kalman filter update. The ODE formulation forms a gradient system with the observations as a cost functional. Besides localization, the new ODE ensemble formulation should also find useful application in the context of nonlinear observation operators and observations that arrive continuously in time.}, language = {en} }