Skip to Content

Theory group software

This page contains brief descriptions of code developed in the group.




The xSPDE code is an extensible Stochastic Partial Differential Equation solver. It is a stochastic toolbox for constructing simulations, which is applicable to many stochastic problems. It has a modular design which can be changed to suit different applications, and includes strategies for calculating errors. At a basic level just one or two lines of input are enough to specify the equation. For advanced users, the entire architecture is open and extensible in numerous ways.

Stochastic equations are equations with random noise terms. They occur in in many fields of science, engineering, economics and other disciplines. xSPDE can solve both ordinary and partial differential stochastic equations. These include partial spatial derivatives, like the Maxwell or Schrödinger equations.

There are many equations of this type, and xSPDE can treat a wide range. It has a configurable functional design. The general structure permits drop-in replacements of the functions provided. Different simulations can be carried out sequentially, to simulate the various stages in an experiment or other process.

The code supports parallelism at both the vector instruction level and at the thread level, using Matlab matrix instructions and the parallel toolbox. It calculates averages of arbitrary functions of any number of complex or real fields. It uses sub-ensemble averaging and extrapolation to obtain accurate error estimates.




PyPi entry:

Reikna is a Python library that contains generalizes and unifies GPGPU code written for various research projects.

It consists of two layers. The low-level layer handles OpenCL/CUDA interface differences, optimal kernel launch parameters, numpy to GPU interfacing (in particular, passing arrays of custom structures to and from GPU). It also includes a temporary memory manager that allows one to greatly reduce the memory footprint of a program if the usage frames of different temporary buffers are known (e.g. it can pack two buffers into one allocation if it knows that they are not used at the same time).

The high-level part is a framework for implementing GPGPU algorithms. The main principle it uses is separating not purely parallel cores of algorithms (that is, the parts that need to exchange data between threads) from purely parallel transformations of the data (like scaling or typecast). It uses the low-level layer functionality to manage temporary buffers needed by algorithms (so that the user does not have to worry about reusing buffers). Currently implemented algorithms are:

  • FFT (a limted version: non power-of-2 sizes are processed much slower) and the related element permutation (fftshift),
  • RNG (counter-based, a port of Random123 algorithm),
  • matrix multiplication,
  • array axis permutation (a generalized version of transposition),
  • reduction with a custom predicate,
  • discrete harmonic transform (an FFT analogue for harmonic trap eigenfunctions).

The main two drawbacks of this library are currently:

  • An algorithm developer has to write some Mako-templated C code, which takes some time to get used to. A transformer from Python functions would be more convenient, but it requires some form of partial evaluation for functions.
  • The API does not really favour heterogeneous or multi-GPU algorithms (one can do those, but they'll have to split and join the data and send it to different threads manually).

Reikna-based ODE/SDE integrator


A Python module based on Reikna that allows one to perform ODE integration on multidimensional fixed grids. Supported integration algorithms are:

The integrator can calculate weak and strong convergence (by propagating with doubled time step, and generating noises accordingly in the case of stochastic equations) and has a simple adaptive step interface (based on the fixed step integration).

Current drawbacks:

  • the user has to write dirft/noise terms in Mako-templated C,
  • the integrator does not split trajectories into subensembles (the sampling error is calculated from the standard deviation over all the trajectories at once), although it can be done manually by a user — there's a function that will merge the results produced by several separate integrator runs and calculate the sampling error from that.




A Python package for numeric simulations of BEC in a harmonic trap. In addition to calculating the dynamics of a condensate, it supports ground state calculation (via imaginary time propagation). The user can collect various predefined observables (population, density, projections and slices of density over given axes, energy) as well as defining custom ones. Both classic GPE and Wigner representation are supported. For the latter, the package can automatically generate noise terms in the equations for given many-body losses.

[an error occurred while processing this directive]