Kassiopeia: A Modern, Extensible C++ Particle Tracking Package

The Kassiopeia particle tracking framework is an object-oriented software package using modern C++ techniques, written originally to meet the needs of the KATRIN collaboration. Kassiopeia features a new algorithmic paradigm for particle tracking simulations which targets experiments containing complex geometries and electromagnetic fields, with high priority put on calculation efficiency, customizability, extensibility, and ease of use for novice programmers. To solve Kassiopeia's target physics problem the software is capable of simulating particle trajectories governed by arbitrarily complex differential equations of motion, continuous physics processes that may in part be modeled as terms perturbing that equation of motion, stochastic processes that occur in flight such as bulk scattering and decay, and stochastic surface processes occuring at interfaces, including transmission and reflection effects. This entire set of computations takes place against the backdrop of a rich geometry package which serves a variety of roles, including initialization of electromagnetic field simulations and the support of state-dependent algorithm-swapping and behavioral changes as a particle's state evolves. Thanks to the very general approach taken by Kassiopeia it can be used by other experiments facing similar challenges when calculating particle trajectories in electromagnetic fields. It is publicly available at https://github.com/KATRIN-Experiment/Kassiopei

A μ-mode BLAS approach for multidimensional tensor-structured problems

In this manuscript, we present a common tensor framework which can be used to generalize one-dimensional numerical tasks to arbitrary dimension d by means of tensor product formulas. This is useful, for example, in the context of multivariate interpolation, multidimensional function approximation using pseudospectral expansions and solution of stiff differential equations on tensor product domains. The key point to obtain an efficient-to-implement BLAS formulation consists in the suitable usage of the mu-mode product (also known as tensor-matrix product or mode-n product) and related operations, such as the Tucker operator. Their MathWorks MATLAB (R)/GNU Octave implementations are discussed in the paper, and collected in the package KronPACK. We present numerical results on experiments up to dimension six from different fields of numerical analysis, which show the effectiveness of the approach

Monolithic multigrid methods for high-order discretizations of time-dependent PDEs

A currently growing interest is seen in developing solvers that couple high-fidelity and higher-order spatial discretization schemes with higher-order time stepping methods for various time-dependent fluid plasma models. These problems are famously known to be stiff, thus only implicit time-stepping schemes with certain stability properties can be used. Of the most powerful choices are the implicit Runge-Kutta methods (IRK). However, they are multi-stage, often producing a very large and nonsymmetric system of equations that needs to be solved at each time step. There have been recent efforts on developing efficient and robust solvers for these systems. We have accomplished this by using a Newton-Krylov-multigrid approach that applies a multigrid preconditioner monolithically, preserving the system couplings, and uses Newton’s method for linearization wherever necessary. We show robustness of our solver on the single-fluid magnetohydrodynamic (MHD) model, along with the (Navier-)Stokes and Maxwell’s equations. For all these, we couple IRK with higher-order (mixed) finiteelement (FEM) spatial discretizations. In the Navier-Stokes problem, we further explore achieving more higher-order approximations by using nonconforming mixed FEM spaces with added penalty terms for stability. While in the Maxwell problem, we focus on the rarely used E-B form, where both electric and magnetic fields are differentiated in time, and overcome the difficulty of using FEM on curved domains by using an elasticity solve on each level in the non-nested hierarchy of meshes in the multigrid method

Computational and Theoretical Developements for (Time Dependent) Density Functional Theory

En esta tesis se presentan avances computacionales y teoricos en la teoria de funcionales de la densidad (DFT) y en la teoria de funcionales de la densidad dependientes del tiempo (TDDFT). Hemos explorado una posible nueva ruta para la mejora de los funcionales de intercambio y correlacion (XCF) en DFT, comprobado y desarrollado propagadores numericos para TDDFT, y aplicado una combinacion de la teoria de control optimo con TDDFT.En los ultimos anos, DFT se ha convertido en el metodo mas utilizado en el area de estructura electronica gracias a su inigualable relacion entre coste y precision. Podemos usar DFT para calcular multitud de propiedades fisicas y quimicas de atomos, moleculas, nanoestructuras, y materia macroscopica. El factor principal que determina la precision que podemos alcanzar usando DFT es el XCF, un objeto desconocido para el cual se han propuesto cientos de aproximaciones distintas. Algunas de estas aproximaciones funcionan correctamente en ciertas situaciones, pero a dia de hoy no existe un XCF que pueda aplicarse con certeza sobre su validez a un sistema arbitrario. Mas aun, no hay una forma sistematica de refinar estos funcionales. Proponemos y exploramos, para sistemas unidimensionales, una nueva manera de estudiarlos y optimizarlos basada en establecer una relacion con la interaccion entre electrones.TDDFT es la extension de DFT a problemas dependientes del tiempo y problemas conestados excitados, y es tambien uno de los metodos mas populares (a veces el unico metodo que se puede poner en practica) en la comunidad de estructura electronica para tratar conellos. De nuevo, la razon detras de su popularidad reside en su relacion precision/coste computacional, que nos permite tratar sistemas mayores y mas complejos. Puede usarse en combinacion con la dinamica de Ehrenfest, un tipo de dinamica molecular no adiabatica.Hemos ido mas alla y hemos combinado TDDFT y la dinamica de Ehrenfest con la teoria de control optimo, creando un instrumento que nos permite, por ejemplo, predecir la forma de los pulsos laser que inducen una explosion de Coulomb en clusters de sodio. A pesar del buen rendimiento computacional de TDDFT en comparacion con otros metodos, hallamos que el coste de estos calculos era bastante elevado.Motivados por este hecho, tambien dedicamos una parte del trabajo de la tesis a la investigacion computacional. En particular, hemos estudiado e implementado familias de propagadores numericos que no se habian examinado en el contexto de TDDFT. Mas concretamente, metodos con varios pasos previos, formulas Runge-Kutta exponenciales, y las expansiones de Magnus sin conmutadores. Finalmente, hemos implementado modificaciones de estas expansiones de Magnus sin conmutadores para la propagacion de las ecuaciones clasico-cuanticas que resultan de la combinacion de la dinamica de Ehrenfest con TDDFT.In this thesis we present computational and theoretical developments for density functional theory (DFT) and time dependent density functional theory (TDDFT). We have explored a new possible route to improve exchange and correlation functionals (XCF) in DFT, tested and developed numerical propagators for TDDFT, and applied a combination of optimal control theory with TDDFT. In recent years, DFT has become the most used method in the electronic structure field thanks to its unparalleled precision/computational cost relationship. We can use DFT to accurately calculate many physical and chemical properties of atoms, molecules, nanostructures, and bulk materials. The main factor that determines the precision that we can obtain using DFT is the XCF, an unknown object for which hundreds of different approximations have been proposed. Some of these approximations work well enough for certain situations, but to this day there is no XCF that can be reliably applied to any arbitrary system. Moreover, there is no clear way for a systematic refinement of these functionals. We propose and explore, for one-dimensional systems, a new way to optimize them, based on establishing a relationship with the electron-electron interaction. TDDFT is the extension of DFT to time-dependent and excited-states problems, and it is also one of the most popular methods (sometimes the only practical one) in the electronic structure community to deal with them. Once again, the reason behind its popularity is its accuracy/computational cost ratio, which allows us to tackle bigger, more complex systems. It can be used in combination with Ehrenfest dynamics, a non-adiabatic type of molecular dynamics. We have furthermore combined both TDDFT and Ehrenfest dynamics with optimal control theory, a scheme that has allowed us, for example, to predict the shapes of the laser pulses that induce a Coulomb explosion in different sodium clusters. Despite the good numerical performance of TDDFT compared to other methods, we found that these computations were still quite expensive. Motivated by this fact, we have also dedicated a part of the thesis work to computational research. In particular, we have studied and implemented families of numerical propagators that had not been tested in the context of TDDFT. More concretely, linear multistep schemes, exponential Runge-Kutta formulas, and commutator-free Magnus expansions. Moreover, we have implemented modifications of these commutator-free Magnus methods for the propagation of the classical-quantum equations that result of combining Ehrenfest dynamics with TDDFT.<br /

Problems in cosmology and numerical relativity

Includes bibliographical references.A generic feature of most inflationary scenarios is the generation of primordial perturbations. Ordinarily, such perturbations can interact with a weak magnetic field in a plasma, resulting in a wide range of phenomena, such as the parametric excitation of plasma waves by gravitational waves. This mechanism has been studied in different contexts in the literature, such as the possibility of indirect detection of gravitational waves through electromagnetic signatures of the interaction. In this work, we consider this concept in the particular case of magnetic field amplification. Specifically, we use non-linear gauge-in variant perturbation theory to study the interaction of a primordial seed magnetic field with density and gravitational wave perturbations in an almost Friedmann-Lemaıtre-Robertson- Walker (FLRW) spacetime with zero spatial curvature. We compare the effects of this coupling under the assumptions of poor conductivity, perfect conductivity and the case where the electric field is sourced via the coupling of velocity perturbations to the seed field in the ideal magnetohydrodynamic (MHD) regime, thus generalizing, improving on and correcting previous results. We solve our equations for long wavelength limits and numerically integrate the resulting equations to generate power spectra for the electromagnetic field variables, showing where the modes cross the horizon. We find that the interaction can seed Electric fields with non-zero curl and that the curl of the electric field dominates the power spectrum on small scales, in agreement with previous arguments. The second focus area of the thesis is the development a stable high order mesh refinement scheme for the solution of hyperbolic partial differential equations. It has now become customary in the field of numerical relativity to couple high order finite difference schemes to mesh refinement algorithms. This approach combines the efficiency of local mesh refinement with the robustness and accuracy of higher order methods. To this end, different modifications of the standard Berger-Oliger adaptive mesh refinement a logarithm have been proposed. In this work we present a new fourth order convergent mesh refinement scheme with sub- cycling in time for numerical relativity applications. One of the distinctive features of our algorithm is that we do not use buffer zones to deal with refinement boundaries, as is currently done in the literature, but explicitly specify boundary data for refined grids instead. We argue that the incompatibility of the standard mesh refinement algorithm with higher order Runge Kutta methods is a manifestation of order reduction phenomena which is caused by inconsistent application of boundary data in the refined grids. Indeed, a peculiar feature of high order explicit Runge Kutta schemes is that they behave like low order schemes when applied to hyperbolic problems with time dependent Dirichlet boundary conditions. We present a new algorithm to deal with this phenomenon and through a series of examples demonstrate fourth order convergence. Our scheme also addresses the problem of spurious reflections that are generated when propagating waves cross mesh refinement boundaries. We introduce a transition zone on refined levels within which the phase velocity of propagating modes is allowed to decelerate in order to smoothly match the phase velocity of coarser grids. We apply the method to test problems involving propagating waves and show a significant reduction in spurious reflections

Exponential integrators: tensor structured problems and applications

The solution of stiff systems of Ordinary Differential Equations (ODEs), that typically arise after spatial discretization of many important evolutionary Partial Differential Equations (PDEs), constitutes a topic of wide interest in numerical analysis. A prominent way to numerically integrate such systems involves using exponential integrators. In general, these kinds of schemes do not require the solution of (non)linear systems but rather the action of the matrix exponential and of some specific exponential-like functions (known in the literature as phi-functions). In this PhD thesis we aim at presenting efficient tensor-based tools to approximate such actions, both from a theoretical and from a practical point of view, when the problem has an underlying Kronecker sum structure. Moreover, we investigate the application of exponential integrators to compute numerical solutions of important equations in various fields, such as plasma physics, mean-field optimal control and computational chemistry. In any case, we provide several numerical examples and we perform extensive simulations, eventually exploiting modern hardware architectures such as multi-core Central Processing Units (CPUs) and Graphic Processing Units (GPUs). The results globally show the effectiveness and the superiority of the different approaches proposed

Programming for Computations - Python: A Gentle Introduction to Numerical Simulations with Python

Numerical simulations; programming; Pytho

Statistical Inference on Dynamical Systems

The ordinary differential equation (ODE) is one representative and popular tool in modeling dynamical systems, which are widely implemented in physics, biology, economics, chemistry and biomedical sciences, etc. Because of the importance of dynamical systems in scientific studies, they are the main focuses of my dissertation. The first chapter of the dissertation is introduction and literature review, which mainly focuses on numerical integration algorithms of ODEs that are difficult to solve analytically, as well as derivative-free optimization algorithms for the so-called inverse problem. The second chapter is on the estimation method based on numerical solvers of differential equations. We start by reviewing the state-of-the-art Gauss-Newton algorithm based method, with the derivation of approximate confidence intervals. Furthermore, we propose and illustrate a method using Differential Evolution along with numerical ODE integration algorithms, as well as a hybrid method to improve the convergence issue for Gauss-Newton algorithm. A numerical comparison study shows the hybrid method is more numerically stable than the traditional Gauss-Newton algorithm based estimation method. In Chapter 3 we propose a novel two-step estimation method based on Fourier basis smoothing and pseudo least square estimator. It is less computationally intensive than methods using numerical ODE integration algorithms, and it works better on periodic or near periodic ODE model functions. In Chapter 4 we expand our study to a population-based hierarchical model to study the correlation between individual features and certain parameter values. Both ML and REML estimation are studied, with more emphasis on REML. An iterative estimation method that incorporates numerical ODE solver into the stochastic approximation EM algorithm for the hierarchical model is proposed and illustrated. Several simulation studies are presented, and a parallel version of the algorithm is implemented as well