Search

A parallel-in-time algorithm based on an augmented Lagrangian approach is proposed to solve four-dimensional variational (4D-Var) data assimilation problems. The assimilation window is divided into multiple sub-intervals that allows to parallelize cost function and gradient computations. Solution continuity equations across interval boundaries are added as constraints. The augmented Lagrangian approach leads to a different formulation of the variational data assimilation problem than weakly constrained 4D-Var. A combination of serial and parallel 4D-Vars to increase performance is also explored. The methodology is illustrated on data assimilation problems with Lorenz-96 and the shallow water models., Comment: 22 Pages

In this article, I will present the postromantic movement, in the hopes that the the aesthetic future of art will belong to those artists who truly value the beauty consecrated for thousands of years by artistic talent. In my opinion, only artistic talent can make the difference between a dilettante and a true professional artist. I long to see the day when artistic value will be recognized, as it used to be, as the result of practical accomplishments not through the possession of a diploma or national ranking of art or resounding critical reviews. I believe that a true artist is not formed only in art schools. A true artist has, deep within his being, an intuitive combination of wish, dream and vision. This gift is only brought to fruition only by those who have the will of saints and the technique of the great masters.

Data assimilation obtains improved estimates of the state of a physical system by combining imperfect model results with sparse and noisy observations of reality. In the four dimensional variational (4D-Var) framework data assimilation is formulated as an optimization problem, which is solved using gradient based optimization methods. The 4D-Var gradient is obtained by forcing the adjoint model with observation increments. The construction of the adjoint model requires considerable development effort. In the continuous approach the adjoint differential equations are discretized. In the discrete approach the numerical solution of the forward equations is differentiated. The two routes lead to different gradients. In this paper we investigate numerically the effect of using discrete and continuous adjoints of the advection equation in chemical transport modeling. Continuous advection adjoints are easily implemented by calling the same advection subroutines as the forward model, with a sign change for the winds, and with rescaling the solution. Discrete advection adjoints involve the differentiation of a nonlinear, monotonic advection scheme like the piecewise parabolic method. The numerical experiments are carried out with CMAQ-ADJ. The results show that, while discrete advection adjoints are more accurate in point-to-point comparisons against finite differences, the continuous adjoints of advection perform better as gradients for optimization in 4D-Var data assimilation.

We present version 3.0 of the atmospheric chemistry box model CAABA/MECCA. In addition to a complete update of the rate coefficients to the most recent recommendations, a number of new features have been added: chemistry in multiple aerosol size bins; automatic multiple simulations reaching steady-state conditions; Monte-Carlo simulations with randomly varied rate coefficients within their experimental uncertainties; calculations along Lagrangian trajectories; mercury chemistry; more detailed isoprene chemistry; tagging of isotopically labeled species. Further changes have been implemented to make the code more user-friendly and to facilitate the analysis of the model results. Like earlier versions, CAABA/MECCA-3.0 is a community model published under the GNU General Public License.

In this paper, the NOx emission scaling factors applied over the 2001 National Emissions Inventory (NEI) are estimated through a four-dimensional variational (4D-Var) approach using SCIAMACHY (Scanning Imaging Absorption spectroMeter for Atmospheric CHartographY) tropospheric NO2 columns measured during summer 2004. In the “top-down” approach, two-month average NO2 columns are assimilated into a regional chemical transport model (CTM), STEM, using different assimilation setups. In a basic setup, NOx emissions are adjusted by assimilating the NO2 columns. A more general setup of emission inversion allows the initial O3 concentrations be adjusted along with the NOx emissions. A final case is set up to assimilate both the NO2 columns and O3 measurement from various platforms while allowing adjustments of both the NOx emissions and the initial O3 concentrations. It is found that the addition of O3 measurements did not improve the NOx emission inversion. With the NOx emission at surface and upper levels being adjusted separately, results from four cases show that the elevated NOx emission reduction ranges from 8.9% to 11.4%, and the surface NOx emission reduction is up to 6.6%. All the cases show NOx emission reduction in Ohio valley and Washington, District of Columbia areas.

Parameter estimation for large systems is a difficult problem, and the solution approaches are computationally expensive. The polynomial chaos approach has been shown to be more efficient than Monte Carlo for quantifying the effects of uncertainties on the system response. This article compares two new computational approaches for parameter estimation based on the polynomial chaos theory for parameter estimation: a Bayesian approach, and an approach using an extended Kalman filter (EKF) to obtain the polynomial chaos representation of the uncertain states and the uncertain parameters. The two methods are applied to a non-linear four-degree-of-freedom roll plane model of a vehicle, in which an uncertain mass with an uncertain position is added on the roll bar. When using appropriate excitations, the results obtained with both approaches are close to the actual values of the parameters, and both approaches can work with noisy measurements. The EKF approach has an advantage over the Bayesian approach: the estimation comes in the form of a posteriori probability densities of the estimated parameters. However, it can yield poor estimations when dealing with non-identifiable systems, and it is recommended to repeat the estimation with different sampling rates in order to verify the coherence of the results with the EKF approach. The Bayesian approach is more robust, can recognize non-identifiability, and use regularization techniques if necessary.

PurposeThe purpose of this paper is to propose a new computational approach for parameter estimation in the Bayesian framework. A posteriori probability density functions are obtained using the polynomial chaos theory for propagating uncertainties through system dynamics. The new method has the advantage of being able to deal with large parametric uncertainties, non‐Gaussian probability densities and nonlinear dynamics.Design/methodology/approachThe maximum likelihood estimates are obtained by minimizing a cost function derived from the Bayesian theorem. Direct stochastic collocation is used as a less computationally expensive alternative to the traditional Galerkin approach to propagate the uncertainties through the system in the polynomial chaos framework.FindingsThe new approach is explained and is applied to very simple mechanical systems in order to illustrate how the Bayesian cost function can be affected by the noise level in the measurements, by undersampling, non‐identifiablily of the system, non‐observability and by excitation signals that are not rich enough. When the system is non‐identifiable and an a priori knowledge of the parameter uncertainties is available, regularization techniques can still yield most likely values among the possible combinations of uncertain parameters resulting in the same time responses than the ones observed.Originality/valueThe polynomial chaos method has been shown to be considerably more efficient than Monte Carlo in the simulation of systems with a small number of uncertain parameters. This is believed to be the first time the polynomial chaos theory has been applied to Bayesian estimation.

This paper discusses the implementation and performance of an array of gas-phase chemistry solvers for the state-of-the-science GEOS-Chem global chemical transport model. The implementation is based on the Kinetic PreProcessor (KPP). Two perl parsers automatically generate the needed interfaces between GEOS-Chem and KPP, and allow access to the chemical simulation code without any additional programming effort. This work illustrates the potential of KPP to positively impact global chemical transport modeling by providing additional functionality as follows. (1) The user can select a highly efficient numerical integration method from an array of solvers available in the KPP library. (2) KPP offers a wide variety of user options for studies that involve changing the chemical mechanism (e.g., a set of additional reactions is automatically translated into efficient code and incorporated into a modified global model). (3) This work provides access to tangent linear, continuous adjoint, and discrete adjoint chemical models, with applications to sensitivity analysis and data assimilation.

On multiprocessors with explicitly managed memory hierarchies (EMM), software has the responsibility of moving data in and out of fast local memories. This task can be complex and error-prone even for expert programmers. Before we can allow compilers to handle this complexity for us, we must identify the abstractions that are general enough to allow us to write applications with reasonable effort, yet specific enough to exploit the vast on-chip memory bandwidth of EMM multi-processors. To this end, we compare two programming models against hand-tuned codes on the STI Cell, paying attention to programmability and performance. The first programming model, Sequoia, abstracts the memory hierarchy as private address spaces, each corresponding to a parallel task. The second, Cellgen, is a new framework which provides OpenMP-like semantics and the abstraction of a shared address space divided into private and shared data. We compare three applications programmed using these models against their hand-optimized counterparts in terms of abstractions, programming complexity, and performance.

In this paper, we discuss the theory of adjoint sensitivity analysis and the method of 4D-Var data assimilation in the context of the Sulfur Transport Eulerian Model (STEM). The STEM atmospheric Chemical Transport Model (CTM) is used to perform adjoint sensitivity analysis and data assimilation over the State of Texas. We first demonstrate the use of adjoint sensitivity analysis for a receptor located at ground level in the Dallas Fort Worth (DFW) area. Simulations are carried out for three 36-h intervals in July 2004. One set of simulations focuses on a passive tracer, and illustrates the influence of meteorological conditions. The other results show the areas of influence associated with DFW ground level ozone, i.e. the areas where changes in precursors (O3, NO2, and HCHO) have the maximum impact on DFW ozone. Next, we employ data assimilation to optimize initial conditions of chemical fields and ground level emissions. We optimize the initial conditions for two episodes on 1 and 16 July 2004. Data assimilation makes use of AirNow ground level observations (for both episodes) and SCHIAMACHY NO2 and HCHO observations from ENVISAT (for the 16 July episode). The re-analyzed chemical fields show considerable improved agreement with AirNow observations for non-initial conditions. We also perform inverse modeling of ground level emissions using data assimilation under an additional smoothness constraint. The results indicate higher NO2 emission levels than in the current emission inventory in the DFW area, and lower emission levels in eastern Texas. The formaldehyde emissions are found to be larger than reported in a localized area near the Gulf Coast, and about at the reported level elsewhere. While the results obtained with the current state-of-the-art tools can help guide tuning of emission inventories, better constraints on the inverse problem are needed to obtain more rigorous quantitative estimates.

We discuss an adaptive resolution system for modeling regional air pollution based on the chemical transport model STEM. The grid adaptivity is implemented using the generic adaptive mesh refinement tool Paramesh, which enables the grid management operations while harnessing the power of parallel computers. The computational algorithm is based on a decomposition of the domain, with the solution in different subdomains being computed with different spatial resolutions. Various refinement criteria that adaptively control the fine grid placement are analyzed to maximize the solution accuracy while maintaining an acceptable computational cost. Numerical experiments in a large-scale parallel setting (~0.5 billion variables) confirm that adaptive resolution, based on a well-chosen refinement criterion, leads to the decrease in spatial error with an acceptable increase in computational time. Fully dynamic grid adaptivity for air quality models is relatively new. We extend previous work on chemical and transport modeling by using dynamically adaptive grid resolution. Advantages and shortcomings of the present approach are also discussed.

This paper discusses several aspects related to the consistency of the discrete adjoints of upwind numerical schemes. Both linear (finite differences, finite volumes) and nonlinear (slope and flux-limited) discretizations of the one-dimensional advection equation are considered. The analysis is focused on uniform meshes and on explicit numerical schemes. We show that the discrete adjoints may lose consistency near the points where upwinding changes, near inflow boundaries where another numerical scheme is employed, and near the locations where the slope/flux limiter is active in the forward simulation. Numerical results are presented to support the theoretical analysis. Copyright © 2007 John Wiley & Sons, Ltd.

Multirate GARK schemes define a multirate extension of GARK schemes, generalized additive Runge-Kutta schemes. In contrast to additive schemes, GARK schemes allow for different stage values as arguments of different components of the right hand side. They introduce additional flexibility when compared to traditional partitioned Runge-Kutta methods, and therefore offer additional opportunities for the development of flexible solvers for systems with multiple scales, or driven by multiple physical processes.

This paper develops an efficient implementation of the ensemble Kalman filter based on a modified Cholesky decomposition for inverse covariance matrix estimation. This implementation is named EnKF-MC. Background errors corresponding to distant model components with respect to some radius of influence are assumed to be conditionally independent. This allows to obtain sparse estimators of the inverse background error covariance matrix. The computational effort of the proposed method is discussed and different formulations based on various matrix identities are provided. Furthermore, an asymptotic proof of convergence with regard to the ensemble size is presented. In order to assess the performance and the accuracy of the proposed method, experiments are performed making use of the Atmospheric General Circulation Model SPEEDY. The results are compared against those obtained using the local ensemble transform Kalman filter (LETKF). Tests are performed for dense observations ($100\%$ and $50\%$ of the model components are observed) as well as for sparse observations (only $12\%$, $6\%$, and $4\%$ of model components are observed). The results reveal that the use of modified Cholesky for inverse covariance matrix estimation can reduce the impact of spurious correlations during the assimilation cycle, i.e., the results of the proposed method are of better quality than those obtained via the LETKF in terms of root mean square error.

Hybrid Monte-Carlo (HMC) sampling smoother is a fully non-Gaussian four-dimensional data assimilation algorithm that works by directly sampling the posterior distribution formulated in the Bayesian framework. The smoother in its original formulation is computationally expensive due to the intrinsic requirement of running the forward and adjoint models repeatedly. Here we present computationally efficient versions of the HMC sampling smoother based on reduced-order approximations of the underlying model dynamics. The schemes developed herein are tested numerically using the shallow-water equations model on Cartesian coordinates. The results reveal that the reduced-order versions of the smoother are capable of accurately capturing the posterior probability density, while being significantly faster than the original full order formulation., 32 pages, 2 figures

An adjoint model for the internationally used Community Multiscale Air Quality (CMAQ) modeling platform of the U.S. EPA is developed. The adjoint version for CMAQ (CMAQ-ADJ) provides the user community with forward (decoupled direct method or DDM) and backward (adjoint) sensitivity analysis capabilities. Current implementation is for gas-phase processes. Discrete adjoints are implemented for all processes with the exception of horizontal advection, for which, because of inherent discontinuities in the advection scheme, the continuous approach is superior. The adjoint of chemistry is constructed by interfacing CMAQ with the kinetic pre-processor, which provides for increased flexibility in the choice of chemical solver and facilitates the implementation of new chemical mechanisms. The adjoint implementation is evaluated both on a process-by-process basis and for the full model. In general, adjoint results show good agreement with brute-force and DDM sensitivities. As expected for a continuous adjoint implementation in a nonlinear scheme, the agreement is not perfect for horizontal transport. Sensitivities of various air quality, public health, and environmental metrics with respect to emissions are calculated using the adjoint method. In order to show applicability to regional climate studies, as an example, the sensitivities of these metrics with respect to local temperatures are calculated.

Data assimilation is the process of integrating observational data and model predictions to obtain an optimal representation of the state of the atmosphere. As more chemical observations in the troposphere are becoming available, chemical data assimilation is expected to play an essential role in air quality forecasting, similar to the role it has in numerical weather prediction (NWP). Considerable progress has been made recently in the development of variational tools for chemical data assimilation. In this paper, we assess the performance of the ensemble Kalman filter (EnKF). Results in an idealized setting show that EnKF is promising for chemical data assimilation.

Data assimilation is the process of integrating observational data and model predictions to obtain an optimal representation of the state of the atmosphere. As more chemical observations in the troposphere are becoming available, chemical data assimilation is expected to play an essential role in air-quality forecasting, similar to the role it has in numerical weather prediction. Considerable progress has been made recently in the development of variational tools for chemical data assimilation. In this paper we assess the performance of the ensemble Kalman filter (EnKF) and compare it with a state-of-the-art 4D-Var approach. We analyse different aspects that affect the assimilation process, and investigate several ways to avoid filter divergence. Results with a real model and real observations show that EnKF is a promising approach for chemical data assimilation. The results also point to several issues on which further research is necessary. Copyright © 2007 Royal Meteorological Society

Data assimilation is the process of integrating observational data and model predictions to obtain an optimal representation of the state of the atmosphere. As more chemical observations in the troposphere are becoming available, chemical data assimilation is expected to play an essential role in air-quality forecasting, similar to the role it has in numerical weather prediction. Considerable progress has been made recently in the development of variational tools for chemical data assimilation. In this paper we implement, and assess the performance of, a localized ‘perturbed-observations’ ensemble Kalman filter (LEnKF). We analyse different settings of the ensemble localization, and investigate the joint assimilation of the state, emissions and boundary conditions. Results with a real model and real observations show that LEnKF is a promising approach for chemical data assimilation. The results also point to several issues on which future research is necessary. Copyright © 2007 Royal Meteorological Society

This work studies reduced order modeling (ROM) approaches to speed up the solution of variational data assimilation problems with large scale nonlinear dynamical models. It is shown that a key requirement for a successful reduced order solution is that reduced order Karush-Kuhn-Tucker conditions accurately represent their full order counterparts. In particular, accurate reduced order approximations are needed for the forward and adjoint dynamical models, as well as for the reduced gradient. New strategies to construct reduced order based are developed for Proper Orthogonal Decomposition (POD) ROM data assimilation using both Galerkin and Petrov-Galerkin projections. For the first time POD, tensorial POD, and discrete empirical interpolation method (DEIM) are employed to develop reduced data assimilation systems for a geophysical flow model, namely, the two dimensional shallow water equations. Numerical experiments confirm the theoretical framework for Galerkin projection. In the case of Petrov-Galerkin projection, stabilization strategies must be considered for the reduced order models. The new reduced order shallow water data assimilation system provides analyses similar to those produced by the full resolution data assimilation system in one tenth of the computational time., Comment: 49 pages, 7 figures

The singular vectors of a chemical transport model are the directions of maximum perturbation growth over a finite time interval. They have proved useful for the estimation of error growth, the initialization of ensemble forecasts, and the optimal placement of adaptive observations. The aim of this paper is to address computational aspects of singular vector analysis for atmospheric chemical transport models. The distinguishing feature of these models is the presence of stiff chemical interactions. A projection approach to preserve the symmetry of the tangent linear–adjoint operator for stiff systems is discussed, and extended to 3D chemical transport simulations. Numerical results are presented for a simulation of atmospheric pollution in East Asia in March 2001. The singular values and the structure of the singular vectors depend on the length of the simulation interval, the meteorological data, the location of the optimization region and the selection of optimization species, the choice of error norms, and the size of the optimization region.

This study explores the use of generalized polynomial chaos theory for modeling complex nonlinear multibody dynamic systems in the presence of parametric and external uncertainty. The polynomial chaos framework has been chosen because it offers an efficient computational approach for the large, nonlinear multibody models of engineering systems of interest, where the number of uncertain parameters is relatively small, while the magnitude of uncertainties can be very large (e.g., vehicle-soil interaction). The proposed methodology allows the quantification of uncertainty distributions in both time and frequency domains, and enables the simulations of multibody systems to produce results with “error bars”. The first part of this study presents the theoretical and computational aspects of the polynomial chaos methodology. Both unconstrained and constrained formulations of multibody dynamics are considered. Direct stochastic collocation is proposed as less expensive alternative to the traditional Galerkin approach. It is established that stochastic collocation is equivalent to a stochastic response surface approach. We show that multi-dimensional basis functions are constructed as tensor products of one-dimensional basis functions and discuss the treatment of polynomial and trigonometric nonlinearities. Parametric uncertainties are modeled by finite-support probability densities. Stochastic forcings are discretized using truncated Karhunen-Loeve expansions. The companion paper “Modeling Multibody Dynamic Systems With Uncertainties. Part II: Numerical Applications” illustrates the use of the proposed methodology on a selected set of test problems. The overall conclusion is that despite its limitations, polynomial chaos is a powerful approach for the simulation of multibody systems with uncertainties.

This study applies generalized polynomial chaos theory to model complex nonlinear multibody dynamic systems operating in the presence of parametric and external uncertainty. Theoretical and computational aspects of this methodology are discussed in the companion paper “Modeling Multibody Dynamic Systems With Uncertainties. Part I: Theoretical and Computational Aspects”.

This paper presents a discretization technique for particle dynamics equation based on piecewise continuous approximations of the solution. The growth terms are re-solved using a Discontinuous Galerkin approach, and the coagulation by a collocation approach. The method fits in the general framework recently proposed by the author. The discontinuous approximation allows better representations of distributions with sharp peaks. Numerical tests reveal that accurate solutions can be obtained with a very small number of size bins. An efficient software package AeroSolve which implements the proposed algorithms was developed.

In this paper we develop the algorithmic tools needed for inverse modeling of aerosol dynamics. Continuous and discrete adjoints of the aerosol dynamic equation are derived, as well as sensitivity coefficients with respect to the coagulation kernel, the growth rate, and the emission and deposition coefficients. Numerical tests performed in the twin experiment framework for a single component model problem show that the initial distributions and the dynamic parameters can be recovered from time series of observations of particle size distributions.

The task of providing an optimal analysis of the state of the atmosphere requires the development of efficient computational tools that facilitate an efficient integration of observational data into models. In a variational approach the data assimilation problem is posed as a minimization problem, which requires the sensitivity (derivatives) of a cost functional with respect to problem parameters. The direct decoupled method has been extensively applied for sensitivity studies of air pollution. Adjoint sensitivity is a complementary approach which efficiently calculates the derivatives of a functional with respect to a large number of parameters. In this paper, we discuss the mathematical foundations of the adjoint sensitivity method applied to air pollution models, and present a complete set of computational tools for performing three-dimensional adjoint sensitivity studies. Numerical examples show that three-dimensional adjoint sensitivity analysis provides information on influence areas, which cannot be obtained solely by an inverse analysis of the meteorological fields. Several illustrative data assimilation results in a twin experiments framework, as well as the assimilation of a real data set are also presented.

When performing dynamic analysis of a constrained mechanical system, a set of index three differential-algebraic equations (DAE) describes the time evolution of the system. The paper presents a state space based method for the numerical solution of the resulting DAE. A subset of so-called independent generalized coordinates, equal in number to the number of degrees of freedom of the mechanical system, is used to express the time evolution of the mechanical system. The second-order state space ordinary differential equations (SSODE) that describe the time variation of independent coordinates are numerically integrated using a Rosenbrock-type formula specialized to second-order systems of differential equations. Rosenbrock methods are known to be efficient for medium accuracy integration of stiff systems; they do not require the solution of non-linear systems for the stage values, and possess optimal linear stability properties for stiff integration. The computation of exact Jacobians needed by Rosenbrock formulas is discussed in the context of multibody systems. The companion paper [19] discusses a choice of method coefficients based on a four-stage L-stable Rosenbrock formula and presents numerical results.

When performing dynamic analysis of a constrained mechanical system, a set of index three differential-algebraic equations (DAE) describes the time evolution of the system. A state-space based method for the numerical solution of the resulting DAE has also been developed. The numerical method uses a linearly implicit time stepping formula of the Rosenbrock type, which is suitable for medium accuracy integration of stiff systems. This paper discusses choices of method coefficients and presents numerical results. For stiff mechanical systems, the proposed algorithm is shown to reduce significantly simulation times when compared to state of the art existent algorithms. The better efficiency is due to the use of an L-stable integrator, and a rigorous and general approach to providing analytical derivatives required by it.

The analysis of comprehensive chemical reactions mechanisms, parameter estimation techniques, and variational chemical data assimilation applications require the development of efficient sensitivity methods for chemical kinetics systems. The new release (KPP-1.2) of the kinetic preprocessor (KPP) contains software tools that facilitate direct and adjoint sensitivity analysis. The direct-decoupled method, built using BDF formulas, has been the method of choice for direct sensitivity studies. In this work, we extend the direct-decoupled approach to Rosenbrock stiff integration methods. The need for Jacobian derivatives prevented Rosenbrock methods to be used extensively in direct sensitivity calculations; however, the new automatic and symbolic differentiation technologies make the computation of these derivatives feasible. The direct-decoupled method is known to be efficient for computing the sensitivities of a large number of output parameters with respect to a small number of input parameters. The adjoint modeling is presented as an efficient tool to evaluate the sensitivity of a scalar response function with respect to the initial conditions and model parameters. In addition, sensitivity with respect to time-dependent model parameters may be obtained through a single backward integration of the adjoint model. KPP software may be used to completely generate the continuous and discrete adjoint models taking full advantage of the sparsity of the chemical mechanism. Flexible direct-decoupled and adjoint sensitivity code implementations are achieved with minimal user intervention. In a companion paper, we present an extensive set of numerical experiments that validate the KPP software tools for several direct/adjoint sensitivity applications, and demonstrate the efficiency of KPP-generated sensitivity code implementations.

The Kinetic PreProcessor KPP was extended to generate the building blocks needed for the direct and adjoint sensitivity analysis of chemical kinetic systems. An overview of the theoretical aspects of sensitivity calculations and a discussion of the KPP software tools is presented in the companion paper. In this work the correctness and efficiency of the KPP generated code for direct and adjoint sensitivity studies are analyzed through an extensive set of numerical experiments. Direct-decoupled Rosenbrock methods are shown to be cost-effective for providing sensitivities at low and medium accuracies. A validation of the discrete–adjoint evaluated gradients is performed against the finite difference estimates. The accuracy of the adjoint gradients is measured using a reference gradient value obtained with a standard direct-decoupled method. The accuracy is studied for both constant step size and variable step size integration of the forward/adjoint model and the consistency between the discrete and continuous adjoint models is analyzed. Applications of the KPP-1.2 software package to direct and adjoint sensitivity studies, variational data assimilation, and parameter identification are considered for the comprehensive chemical mechanism SAPRC-99.

This paper deals with the numerical solution of systems of differential equations with a stiff part and a non-stiff one, typically arising from the semi-discretization of certain partial differential equations models. It is illustrated the construction and analysis of highly stable and high-stage order implicit-explicit (IMEX) methods based on diagonally implicit multistage integration methods (DIMSIMs), a subclass of general linear methods (GLMs). Some examples of methods with optimal stability properties are given. Finally numerical experiments confirm the theoretical expectations.

With organic photovoltaic (OPV) technology moving towards commercialisation, high-throughput analytical techniques are required to study the nanoscale morphology of OPV blends. We demonstrate a low-voltage backscattered electron imaging technique in the SEM that combines a solid-state backscattered electron detector with stage biasing to produce a rapid overview of the phase-separated surface morphology of an organic photovoltaic (P3HT:PCBM) blend. Aspects of obtaining the best possible results from the technique are discussed along with the possibility of probing the sub-surface morphology by altering the primary electron beam landing energy.

Implicit-explicit (IMEX) time stepping methods can efficiently solve differential equations with both stiff and nonstiff components. We have recently proposed new implicit-explicit methods of general linear type (IMEX-GLMs) in [11, 10]. It was shown that no additional coupling order conditions are needed, and consequently that GLMs offer an excellent framework for the construction of multi-method integration algorithms. In this paper we investigate the application of IMEX GLMs to solve reaction-diffusion problems. Numerical results show that IMEX GLMs are competitive for these applications.

The kinetic preprocessor (KPP) is a software tool that assists the computer simulation of chemical kinetic systems. The concentrations of a chemical system evolve in time according to the differential law of mass action kinetics. A computer simulation requires the implementation of the differential system and its numerical integration in time. KPP translates a specification of the chemical mechanism into fortran or c simulation code that implement the concentration time derivative function and its Jacobian, together with a suitable numerical integration scheme. Sparsity in Jacobian is carefully exploited in order to obtain computational efficiency. KPP incorporates a library with several widely used atmospheric chemistry mechanisms and users can add their own chemical mechanisms to the library. KPP also includes a comprehensive suite of stiff numerical integrators. The KPP development environment is designed in a modular fashion and allows for rapid prototyping of new chemical kinetic schemes as well as new numerical integration methods.

PAQMSG is an MPI-based, Fortran 90 communication library for the parallelization of air quality models (AQMs) on structured grids. It consists of distribution, gathering and repartitioning routines for different domain decompositions implementing a master–worker strategy. The library is architecture and application independent and includes optimization strategies for different architectures. This paper presents the library from a user perspective. Results are shown from the parallelization of STEM-III on Beowulf clusters. The PAQMSG library is available on the web. The communication routines are easy to use, and should allow for an immediate parallelization of existing AQMs. PAQMSG can also be used for constructing new models.

This paper presents a direct approach to solving the aerosol coagulation equation. Newton–Cotes formulas are used to discretize the integral terms, and the semi-discrete system is built using collocation. A semi-implicit Gauss–Seidel time integration method is employed. The approach generalizes the semi-implicit method of Jacobson and is more accurate.

Implicit---explicit (IMEX) time stepping methods can efficiently solve differential equations with both stiff and nonstiff components. IMEX Runge---Kutta methods and IMEX linear multistep methods have been studied in the literature. In this paper we study new implicit---explicit methods of general linear type. We develop an order conditions theory for high stage order partitioned general linear methods (GLMs) that share the same abscissae, and show that no additional coupling order conditions are needed. Consequently, GLMs offer an excellent framework for the construction of multi-method integration algorithms. Next, we propose a family of IMEX schemes based on diagonally-implicit multi-stage integration methods and construct practical schemes of order up to three. Numerical results confirm the theoretical findings.

Sensitivity analysis of multibody systems is essential for several applications, such as dynamics-based design optimization. Dynamic sensitivities, when needed, are often calculated by means of finite differences. This procedure is computationally expensive when the number of parameters is large, and numerical errors can severely limit its accuracy. This paper explores several analytical approaches to perform sensitivity analysis of multibody systems. Direct and adjoint sensitivity equations are developed in the context of Maggi's formulation of multibody dynamics equations. The approach can be generalized to other formulations of multibody dynamics as systems of ordinary differential equations (ODEs). The sensitivity equations are validated numerically against the third party code fatode and against finite difference solutions with real and complex perturbations.

Linearly implicit Runge-Kutta methods with approximate matrix factorization can solve efficiently large systems of differential equations that have a stiff linear part, e.g. reaction-diffusion systems. However, the use of approximate factorization usually leads to loss of accuracy, which makes it attractive only for low order time integration schemes. This paper discusses the application of approximate matrix factorization with high order methods; an inexpensive correction procedure applied to each stage allows to retain the high order of the underlying linearly implicit Runge-Kutta scheme. The accuracy and stability of the methods are studied. Numerical experiments on reaction-diffusion type problems of different sizes and with different degrees of stiffness illustrate the efficiency of the proposed approach.

This work presents a novel motion planning framework, rooted in nonlinear programming theory, that treats uncertain fully and underactuated dynamical systems described by ordinary differential equations. Uncertainty in multibody dynamical systems comes from various sources, such as system parameters, initial conditions, sensor and actuator noise, and external forcing. Treatment of uncertainty in design is of paramount practical importance because all real-life systems are affected by it, and poor robustness and suboptimal performance result if it is not accounted for in a given design. In this work uncertainties are modeled using generalized polynomial chaos and are solved quantitatively using a least-square collocation method. The computational efficiency of this approach enables the inclusion of uncertainty statistics in the nonlinear programming optimization process. As such, the proposed framework allows the user to pose, and answer, new design questions related to uncertain dynamical systems. Specifically, the new framework is explained in the context of forward, inverse, and hybrid dynamics formulations. The forward dynamics formulation, applicable to both fully and underactuated systems, prescribes deterministic actuator inputs that yield uncertain state trajectories. The inverse dynamics formulation is the dual to that of forward dynamics, and is only applicable to fully actuated systems; deterministic state trajectories are prescribed and yield uncertain actuator inputs. The inverse dynamics formulation is more computationally efficient as it requires only algebraic evaluations and completely avoids numerical integration. Finally, the hybrid dynamics formulation is applicable to underactuated systems where it leverages the benefits of inverse dynamics for actuated joints and forward dynamics for unactuated joints; it prescribes actuated state and unactuated input trajectories that yield uncertain unactuated states and uncertain actuated inputs. The benefits of the ability to quantify uncertainty when planning the motion of multibody dynamic systems are illustrated through several case studies. The resulting designs determine optimal motion plans—subject to deterministic and statistical constraints—for all possible systems within the probability space.

Chemical kinetics conserves mass and renders nonnegative solutions; a good numerical simulation would ideally produce mass-balanced, positive concentration vectors. Many time-stepping methods are mass conservative; however, unconditional positivity restricts the order of a traditional method to one. The projection method presented in this paper ensures mass conservation and positivity. First, a numerical approximation is computed with one step of a mass-preserving traditional scheme. If there are negative components, the nearest vector in the reaction simplex is found by solving a quadratic optimization problem; this vector is shown to better approximate the true solution. A simpler version involves just one projection step and stabilizes the reaction simplex. This technique works best when the underlying time-stepping scheme favors positivity. Projected methods are more accurate than clipping and allow larger time steps for kinetic systems which are unstable outside the positive quadrant.

In the past decade the variational method has been successfully applied in data assimilation problems for atmospheric chemistry models. In 4D-var data assimilation, a minimization algorithm is used to find the set of control variables which minimizes the weighted least squares distance between model predictions and observations over the assimilation window. Using the adjoint method, the gradient of the cost function can be computed fast, at the expense of few function evaluations, making the optimization process very efficient. For large-scale models, the high storage requirements and the difficulty of implementing the adjoint code when sophisticated integrators are used to solve the stiff chemistry make the assimilation a very intensive computational process. If the sparse structure of the chemical models is carefully exploited, Rosenbrock methods have been proved to be reliable chemistry solvers because of their outstanding stability properties and conservation of the linear invariants of the system. In this paper we present an efficient implementation of the adjoint code for the Rosenbrock type methods, which can reduce the storage requirements of the forward model and is suitable for automatization. The adjoint code is completely generated using symbolic preprocessing and automatic differentiation tools which allow flexibility and require minimal user intervention.

Discusses computational challenges in air quality modelling (as viewed by the authors). The focus of the paper will be on Di, the “current” state‐of‐affairs. Owing to limitation of space the discussion will focus on only a few aspects of air quality modelling: i.e. chemical integration, sensitivity analysis and computational framework, with particular emphasis on aerosol issues.

Numerical resonance artifacts have become recognized recently as a limiting factor to increasing the timestep in multiple-timestep (MTS) biomolecular dynamics simulations. At certain timesteps correlated to internal motions (e.g., 5 fs, around half the period of the fastest bond stretch,Tmin), visible inaccuracies or instabilities can occur. Impulse-MTS schemes are vulnerable to these resonance errors since large energy pulses are introduced to the governing dynamics equations when the slow forces are evaluated. We recently showed that such resonance artifacts can be masked significantly by applying extrapolative splitting to stochastic dynamics. Theoretical and numerical analyses of force-splitting integrators based on the Verlet discretization are reported here for linear models to explain these observations and to suggest how to construct effective integrators for biomolecular dynamics that balance stability with accuracy. Analyses for Newtonian dynamics demonstrate the severe resonance patterns of the Impulse splitting, with this severity worsening with the outer timestep, ?t; Constant Extrapolation is generally unstable, but the disturbances do not grow with ?t. Thus, the stochastic extrapolative combination can counteract generic instabilities and largely alleviate resonances with a sufficiently strong Langevin heat-bath coupling (?), estimates for which are derived here based on the fastest and slowest motion periods. These resonance results generally hold for nonlinear test systems: a water tetramer and solvated protein. Proposed related approaches such as Extrapolation/Correction and Midpoint Extrapolation work better than Constant Extrapolation only for timesteps less thanTmin/2. An effective extrapolative stochastic approach for biomolecules that balances long-timestep stability with good accuracy for the fast subsystem is then applied to a biomolecule using a three-class partitioning: the medium forces are treated byMidpoint Extrapolationvia position Verlet, and the slow forces are incorporated byConstant Extrapolation. The resulting algorithm (LN) performs well on a solvated protein system in terms of thermodynamic properties and yields an order of magnitude speedup with respect to single-timestep Langevin trajectories. Computed spectral density functions also show how the Newtonian modes can be approximated by using a small ? in the range of 5?20 ps?1.

Purpose: To quantify the dosimetric effects of patient positioning uncertainties on PTV coverage and normal tissue sparing and find the tolerance for IMRT and SBRT positioning errors.Method and Materials: Patient positioning uncertainties were represented by isocenter shifts. Three head and neck IMRT, three prostate IMRT, and two lung SBRT plans were included in this study. The potential positioning uncertainties were introduced by shifting the isocenters by 0.1 cm, 0.3 cm, 0.5 cm, .0 cm, and 1.5 cm in Eclipse. Actual isocenter shifts (clinical shift) based on OBI/CBCT were also applied to the plans. The plans with isocenter shifts were recalculated in Eclipse and the results were compared with the original plans. Results: The dose errors increased with isocenter shifts. For head and neck patients, the PTV D95 dropped by 3.3% and the cord dose increased by 3.0% for 3-mm shift, but D95 dropped 6.8 % with 5 mm shift, and spinal cord dose increased by 6.9%. For prostate, 3-mm shift reduced 3.4% in D95 and increased 1.2% rectum dose and 2.7% bladder dose, while 5 mm shift caused 9.8% of reduction in D95 and 14.2% increase in bladder dose, but 7.1% decrease in rectum mean dose. For lung SBRT, 1-mm shift caused 2.2% decrease in D95 and 0.7% decrease in spinal cord, but 3 mm shifts led to 14.9% of decrease in D95, while spinal cord dose also decreased by 1.5%.ConclusionThe patient positioning uncertainties reduced the PTV coverage and changed the dose to the normal tissues. However, the significance of dosimetric effect of patient positioning uncertainties varied with different tumor sites. This study suggested that the positioning uncertainty should be within 3 mm for regular IMRT treatment, and 1 mm for SBRT to avoid significant dose errors.

This paper develops a new class of exponential-type integrators where all the matrix exponentiations are performed in a single Krylov space of low dimension. The new family, called Lightly Implicit Krylov-Exponential (LIKE), is well suited for solving large scale systems of ODEs or semi-discrete PDEs. The time discretization and the Krylov space approximation are treated as a single computational process, and the Krylov space properties are an integral part of the new LIKE order condition theory developed herein. Consequently, LIKE methods require a small number of basis vectors determined solely by the temporal order of accuracy. The subspace size is independent of the ODE under consideration, and there is no need to monitor the errors in linear system solutions at each stage. Numerical results illustrate the favorable properties of new family of methods.

The optimization of the dynamic response of multibody dynamic systems is a complex and open problem. It relies on using the equations of motion of the system, which amplifies the level of complexity of the problem substantially, compared to other types of optimization. In the context of this kind of optimization, the sensitivity analysis of the dynamic response of the system is a key element. Two main techniques are currently available for the sensitivity analysis of the response of a dynamical system: the direct differentiation method and the adjoint variable method. In this work, different formulations of the equations of motion with dependent coordinates are employed and their sensitivity equations obtained. Direct and adjoint sensitivities are convenient for different types of problems but both methods require accurate derivatives of the equations of motion considered. Both approaches are employed in this work; the direct and adjoint sensitivity equations are obtained for index-3 differential-algebraic equations (DAE), index-1 DAE, and penalty formulations. The adjoint sensitivity for the penalty formulations introduced here is completely novel.

This paper develops a new class of Rosenbrock-type integrators based on a Krylov space solution of the linear systems. The new family, called Rosenbrock-Krylov (Rosenbrock-K), is well suited for solving large scale systems of ODEs or semi-discrete PDEs. The time discretization and the Krylov space approximation are treated as a single computational process, and the Krylov space properties are an integral part of the new Rosenbrock-K order condition theory developed herein. Consequently, Rosenbrock-K methods require a small number of basis vectors determined solely by the temporal order of accuracy. The subspace size is independent of the ODE under consideration, and there is no need to monitor the errors in linear system solutions at each stage. Numerical results show favorable properties of Rosenbrock-K methods when compared to current Rosenbrock and Rosenbrock-W schemes.