We have defined and analyzed a semi-Toeplitz preconditioner for time-dependent and steady-state convection-diffusion problems. Analytic expressions for the eigenvalues of the preconditioned systems are obtained. An asymptotic analysis shows that the eigenvalue spectrum of the time-dependent problem is reduced to two eigenvalues when the number of grid points go to infinity. The numerical experiments sustain the results of the theoretical analysis, and the preconditioner exhibits a robust behavior for stretched grids. A semi-Toeplitz preconditioner for the linearized Navier-Stokes equations for compressible flow is proposed and tested. The preconditioner is applied to the linear system of equations to be solved in each time step of an implicit method. The equations are solved with flat plate boundary conditions and are linearized around the Blasius solution. The grids are stretched in the normal direction to the plate and the quotient between the time step and the space step is varied. The preconditioner works well in all tested cases and outperforms the method without preconditioning both in number of iterations and execution time.

New preconditioning techniques for the iterative solution of systems of equations arising from discretizations of partial differential equations are considered. Fundamental solutions, both of differential and difference operators, are used as kernels in discrete, truncated convolution operators. The intention is to approximate inverses of difference operators that arise when discretizing the differential equations. The approximate inverses are used as preconditioners. The technique using fundamental solutions of differential operators is applied to scalar problems in two dimensions, and grid independent convergence is obtained for a first order differential equation. The problem of computing fundamental solutions of difference operators is considered, and we propose a new algorithm. It can be used also when the symbol of the difference operator is not invertible everywhere, and it is applicable in two or more dimensions. Fundamental solutions of difference operators are used to construct preconditioners for non-linear systems of difference equations in two dimensions. Grid independent convergence is observed for two standard finite difference discretizations of the Euler equations in a non-axisymmetric duct.

Two unstructured finite volume solvers for the Maxwell equations in 2D and 3D are introduced. The solvers are a generalization of FD–TD to unstructured grids and they use a third-order staggered Adams–Bashforth scheme for time discretization. Analysis and experiments of this time integrator reveal that we achieve a long term stable solution on general triangular grids. A Fourier analysis shows that the 2D solver has excellent dispersion characteristics on uniform triangular grids. In 3D a spatial filter of Laplace type is introduced to enable long simulations without suffering from late time instability. The recursive convolution method proposed by Luebbers et al. to extend FD–TD to permit frequency dispersive materials is here generalized to the 3D solver. A better modelling of materials which have a strong frequency dependence in their constitutive parameters is obtained through the use of a general material model. The finite volume solvers are not intended to be stand-alone solvers but one part in two hybrid solvers with FD–TD. The numerical examples in 2D and 3D demonstrate that the hybrid solvers are superior to stand-alone FD–TD in terms of accuracy and efficiency.

We have developed an iterative solver for the Moment Method. It computes a matrix–vector product with the multilevel Fast Multipole Method, which makes the method scale with the number of unknowns. The iterative solver is of Block Quasi-Minimum Residual type and can handle several right-hand sides at once. The linear system is preconditioned with a Sparse Approximate Inverse, which is modified to handle dense matrices. The solver is parallelized on shared memory machines using OpenMP. To verify the method some tests are conducted on varying geometries. We use simple geometries to show that the method works. We show that the method scales on several processors of a shared memory machine. To prove that the method works for real life problems, we do some tests on large scale aircrafts. The largest test is a one million unknown simulation on a full scale model of a fighter aircraft.

The time-dependent Schrödinger equation models the quantum nature of molecular processes. Numerical simulations of these models help in understanding and predicting the outcome of chemical reactions. In this thesis, several numerical algorithms for evolving the Schrödinger equation with an explicitly time-dependent Hamiltonian are studied and their performance is compared for the example of a pump-probe and an interference experiment for the rubidium diatom. For the important application of interaction dynamics between a molecule and a time-dependent field, an efficient fourth order Magnus-Lanczos propagator is derived. Error growth in the equation is analyzed by means of a posteriori error estimation theory and the self-adjointness of the Hamiltonian is exploited to yield a low-cost global error estimate for numerical time evolution. Based on this theory, an h,p-adaptive Magnus-Lanczos propagator is developed that is capable to control the global error. Numerical experiments for various model systems (including a three dimensional model and a dissociative configuration) show that the error estimate is effective and the number of time steps needed to meet a certain accuracy is reduced due to adaptivity. Moreover, the thesis proposes an efficient numerical optimization framework for the design of femtosecond laser pulses with the aim of manipulating chemical reactions. This task can be formulated as an optimal control problem with the electric field of the laser being the control variable. In the algorithm described here, the electric field is Fourier transformed and it is optimized over the Fourier coefficients. Then, the frequency band is narrowed which facilitates the application of a quasi-Newton method. Furthermore, the restrictions on the frequency band make sure that the optimized pulse can be realized by the experimental equipment. A numerical comparison shows that the new method can outperform the Krotov method, which is a standard scheme in this field.

This thesis discusses the numerical approximation of flow problems, in particular the large eddy simulation of turbulent flow and the simulation of laminar immiscible two-phase flow. The computations for both applications are performed with a coupled solution approach of the Navier-Stokes equations discretized with the finite element method. Firstly, a new implementation strategy for large eddy simulation of turbulent flow is discussed. The approach is based on the variational multiscale method, where scale ranges are separated by variational projection. The method uses a standard Navier-Stokes model for representing the coarser of the resolved scales, and adds a subgrid viscosity model to the smaller of the resolved scales. The scale separation within the space of resolved scales is implemented in a purely algebraic way based on a plain aggregation algebraic multigrid restriction operator. A Fourier analysis underlines the importance of projective scale separations and that the proposed model does not affect consistency of the numerical scheme. Numerical examples show that the method provides better results than other state-of-the-art methods for computations at low resolutions. Secondly, a method for modeling laminar two-phase flow problems in the vicinity of contact lines is proposed. This formulation combines the advantages of a level set model and of a phase field model: Motion of contact lines and imposition of contact angles are handled like for a phase field method, but the computational costs are similar to the ones of a level set implementation. The model is realized by formulating the Cahn-Hilliard equation as an extension of a level set model. The phase-field specific terms are only active in the vicinity of contact lines. Moreover, various aspects of a conservative level set method discretized with finite elements regarding the accuracy of force balance and prediction in jump of pressure between the inside and outside of a circular bubble are tested systematically. It is observed that the errors in velocity are mostly due to inaccuracies in curvature evaluation, whereas the errors in pressure prediction mainly come from the finite width of the interface. The error in both velocity and pressure decreases with increasing number of mesh points.

Chemical dissociation processes are important in quantum dynamics. Such processes can be investigated theoretically and numerically through the time-dependent Schrödinger equation, which gives a quantum mechanical description of molecular dynamics. This thesis discusses the numerical simulation of chemical reactions involving dissociation. In particular, an accurate boundary treatment in terms of artificial, absorbing boundaries of the computational domain is considered. The approach taken here is based on the perfectly matched layer technique in a finite difference framework. The errors introduced due to the perfectly matched layer can be divided into two categories, the modeling error from the continuous model and numerical reflections that arise for the discretized problem. We analyze the different types of errors using plane wave analysis, and parameters of the perfectly matched layer are optimized so that the modeling error and the numerical reflections are of the same order. The level of accuracy is determined by estimating the order of the spatial error in the interior domain. Numerical calculations show that this procedure enables efficient calculations within a given accuracy. We apply our perfectly matched layer to a three-state system describing a one-dimensional IBr molecule subjected to a laser field and to a two-dimensional model problem treating dissociative adsorbtion and associative desorption of an H2 molecule on a solid surface. Comparisons made to standard absorbing layers in chemical physics prove our approach to be efficient, especially when high accuracy is of importance. / eSSENCE

Numerical simulation of propagating waves in unbounded spatial domains is a challenge common to many branches of engineering and applied mathematics. Perfectly matched layers (PML) are a novel technique for simulating the absorption of waves in open domains. The equations modeling the dynamics of phenomena of interest are usually posed as differential equations (or integral equations) which must be solved at every time instant. In many application areas like general relativity, seismology and acoustics, the underlying equations are systems of second order hyperbolic partial differential equations. In numerical treatment of such problems, the equations are often rewritten as first order systems and are solved in this form. For this reason, many existing PML models have been developed for first order systems. In several studies, it has been reported that there are drawbacks with rewriting second order systems into first order systems before numerical solutions are obtained. While the theory and numerical methods for first order systems are well developed, numerical techniques to solve second order hyperbolic systems is an on-going research. In the first part of this thesis, we construct PML equations for systems of second order hyperbolic partial differential equations in two space dimensions, focusing on the equations of linear elasto-dynamics. One advantage of this approach is that we can choose auxiliary variables such that the PML is strongly hyperbolic, thus strongly well-posed. The second is that it requires less auxiliary variables as compared to existing first order formulations. However, in continuum the stability of both first order and second order formulations are linearly equivalent. A turning point is in numerical approximations. We have found that if the so-called geometric stability condition is violated, approximating the first order PML with standard central differences leads to a high frequency instability for any given resolution. The second order discretization behaves much more stably. In the second order setting instability occurs only if unstable modes are well resolved. The second part of this thesis discusses the construction of PML equations for the time-dependent Schrödinger equation. From mathematical perspective, the Schrödinger equation is unique, in the sense that it is only first order in time but second order in space. However, with slight modifications, we carry over our ideas from the hyperbolic systems to the Schrödinger equations and derive a set of asymptotically stable PML equations. The new model can be viewed as a modified complex absorbing potential (CAP). The PML model can easily be adapted to existing codes developed for CAP by accurately discretizing the auxiliary variables and appending them accordingly. Numerical experiments are presented illustrating the accuracy and absorption properties of the new PML model. We are hopeful that the results obtained in this thesis will find useful applications in time-dependent wave scattering calculations.

The aim of the project is to construct, analyse and implement fast and reliable numerical solution methods to simulate multi-phase flow, modeled by a coupled system consisting of the time-dependent Cahn-Hilliard and incompressible Navier-Stokes equations with variable viscosity and variable density. This thesis mainly discusses the efficient solution methods for the latter equations aiming at constructing preconditioners, which are numerically and computationally efficient, and robust with respect to various problem, discretization and method parameters. In this work we start by considering the stationary Navier-Stokes problem with constant viscosity. The system matrix arising from the finite element discretization of the linearized Navier-Stokes problem is nonsymmetric of saddle point form, and solving systems with it is the inner kernel of the simulations of numerous physical processes, modeled by the Navier-Stokes equations. Aiming at reducing the simulation time, in this thesis we consider iterative solution methods with efficient preconditioners. When discretized with the finite element method, both the Cahn-Hilliard equations and the stationary Navier-Stokes equations with constant viscosity give raise to linear algebraic systems with nonsymmetric matrices of two-by-two block form. In Paper I we study both problems and apply a common general framework to construct a preconditioner, based on the matrix structure. As a part of the general framework, we use the so-called element-by-element Schur complement approximation. The implementation of this approximation is rather cheap. However, the numerical experiments, provided in the paper, show that the preconditioner is not fully robust with respect to the problem and discretization parameters, in this case the viscosity and the mesh size. On the other hand, for not very convection-dominated flows, i.e., when the viscosity is not very small, this approximation does not depend on the mesh size and works efficiently. Considering the stationary Navier-Stokes equations with constant viscosity, aiming at finding a preconditioner which is fully robust to the problem and discretization parameters, in Paper II we turn to the so-called augmented Lagrangian (AL) approach, where the linear system is transformed into an equivalent one and then the transformed system is iteratively solved with the AL type preconditioner. The analysis in Paper II focuses on two issues, (1) the influence of a scalar method parameter (a stabilization constant in the AL method) on the convergence rate of the preconditioned method and (2) the choice of a matrix parameter for the AL method, which involves an approximation of the inverse of the finite element mass matrix. In Paper III we consider the stationary Navier-Stokes problem with variable viscosity. We show that the known efficient preconditioning techniques in particular, those for the AL method, derived for constant viscosity, can be straightforwardly applicable also in this case. One often used technique to solve the incompressible Navier-Stokes problem with variable density is via operator splitting, i.e., decoupling of the solutions for density, velocity and pressure. The operator splitting technique introduces an additional error, namely the splitting error, which should be also considered, together with discretization errors in space and time. Insuring the accuracy of the splitting scheme usually induces additional constrains on the size of the time-step. Aiming at fast numerical simulations and using large time-steps may require to use higher order time-discretization methods. The latter issue and its impact on the preconditioned iterative solution methods for the arising linear systems are envisioned as possible directions for future research. When modeling multi-phase flows, the Navier-Stokes equations should be considered in their full complexity, namely, the time-dependence, variable viscosity and variable density formulation. Up to the knowledge of the author, there are not many studies considering all aspects simultaneously. Issues on this topic, in particular on the construction of efficient preconditioners of the arising matrices need to be further studied.

An adjoint system of the Euler equations of gas dynamics is derived in order to solve aerodynamic shape optimization problems with gradient-based methods. The derivation is based on the fully discrete flow model and involves differentiation and transposition of the system of equations obtained by an unstructured and node-centered finite-volume discretization. Solving the adjoint equations allows an efficient calculation of gradients, also when the subject of optimization is described by hundreds or thousands of design parameters. Such a fine geometry description may cause wavy or otherwise irregular designs during the optimization process. Using the one-to-one mapping defined by a Poisson problem is a known technique that produces smooth design updates while keeping a fine resolution of the geometry. This technique is extended here to combine the smoothing effect with constraints on the geometry, by defining the design updates as solutions of a quadratic programming problem associated with the Poisson problem. These methods are applied to airfoil shape optimization for reduction of the wave drag, that is, the drag caused by gas dynamic effects that occur close to the speed of sound. A second application concerns airfoil design optimization to delay the laminar-to-turbulent transition point in the boundary layer in order to reduce the drag. The latter application has been performed by the author with collaborators, also using gradient-based optimization. Here, the growth of convectively unstable disturbances are modeled by successively solving the Euler equations, the boundary layer equations, and the parabolized stability equations.