Home

The cuPSS library is in a very early stage of development but is already capable of solving a wide range of coupled stochastic PDEs.

It does so by numerically integrating PDEs in Fourier space, and thus is particularly good for problems with high order derivatives in which real space formulations might lead to instabilities due to the discretization of gradients.

cuPSS runs natively in CUDA, where 1D or 2D problems are treated as textures. This makes it potentially faster than other packages without this type of hardware acceleration.

Installation

Dependencies

cuPSS relies on FFTW3 and several of the standard CUDA toolkit libraries. Install FFTW3 with single floating point precision. Download from fftw3's website here. Decompress in a directory and run

./configure --enable-float
make
sudo make install

To install the CUDA toolkit, follow the instructions in their website

Make sure the toolkit installed matches the version of the NVIDIA drivers. CUDA will work and run sometimes with mismatching versions, without making it clear there are any errors. Run nvidia-smi to check the CUDA Version, on the top right, and download the matching CUDA toolkit.

A piece of cuPSS relies on the standard filesystem header file to create data directories. This is standard in c++17 and beyond, so compile cuPSS solvers under this standard.

Quickstart

Once the dependencies have been installed, cuPSS can be tried by cloning the repository

git clone https://github.com/fcaballerop/cuPSS.git

Compiling cuPSS, installing system-wide

To compile and install system-wide, CMake can be used, run the following commands in the cuPSS directory.

mkdir build
cd build
cmake ../
cmake --build .
sudo cmake --install .
cd ..

Try one of the example solvers, compile linking the appropriate libraries

nvcc examples/cahn-hilliard.cpp -lcufft -lfftw3f -lcurand -lcupss -O2 -o cahn-hilliard
./cahn-hilliard

The output files will be in a directory called data, in the directory where the executable is. These files can be analysed or plotted with any software.

Local compilation

If CMake is not available, or system-wide installation is not possible, cuPSS can be compiled locally int a static library to then be linked to any solver. To do so, the source files can be compiled to object files and packed into a library by running

cd src
nvcc -c *cu *cpp -O2
ar rcs libcupss.a *.o
cd ..

The library will be in the directory src/, and can be linked to any solver. The header file to include will then be in the inc directory. Try compiling an example solver

nvcc examples/cahn-hilliard.cpp -Lsrc/ -lcufft -lfftw3f -lcurand -lcupss -O2 -o cahn-hilliard
./cahn-hilliard

The executable will run and output data in directory called 'data' located where the executable is.

What can it solve?

cuPSS can solve a system of $N$ coupled first order PDEs. If these $N$ PDEs describe the evolution of of $N$ fields $\lbrace\phi_i\rbrace_{i=1}^N$, then the system of equations can be written as

$$ \begin{aligned} \partial_t \phi_1(k,t) &= F_1[\lbrace\phi_i\rbrace,\lbrace\psi_i\rbrace] + \eta_1 \\ \vdots & \\ \partial_t \phi_N(k,t) &= F_N[\lbrace\phi_i\rbrace,\lbrace\psi_i\rbrace] + \eta_N \end{aligned} $$

where each $F_i$ represents the time evolution of each field $\phi_i$, and we allow for the presence of other nondynamic fields $\psi_i$. Likewise, the equations can be deterministic or stochastic, with noise terms $\eta_i$. As a spectral solver, it integrates the equations of motion in Fourier space and thus the equations are written for the fields in Fourier space ( $\phi(k,t)$ ).

Nonlinearities are calculated pseudospectrally to avoid explicitly calculating convolutions in Fourier space.

See tutorials for example on how to write solvers using cuPSS.

How does it solve equations?

The equations are integrated in time by simple Euler steps in Fourier space. Linear terms can be integrated through forward or backward steps. Nonlinear terms are integrated through forward steps. Nonlinear terms are integrated pseudospectrally, meaning, each field that is multiplied in the nonlinear term will be transformed back to real space, then the nonlinearity is computed by multiplying its components in real space, and then the result is transformed back to real space. Before the first transformation to real space, fields are dealiased. Explicitly, a nonlinear term of the form $\partial_t \phi = \ldots + L \prod_i \phi_i$, with some operator $L$ applied to the whole term is calculated in Fourier space as:

$$ \frac{\phi_{t+\Delta t}(k) - \phi_t(k)}{\Delta t} = \ldots + \tilde L \mathcal F\left[\prod_i \mathcal F^{-1}[\phi_i]\right] $$

where $\mathcal F$ and $\mathcal F^{-1}$ are forward and backward Fourier transforms respectively. For instance, a second order nonlinearity of the form $\partial_t\phi(r,t) = \ldots + \nabla^2\left((\nabla\phi)^2\right)$ would be calculated as

$$ \frac{\phi_{t+\Delta t}(k) - \phi_t(k)}{\Delta t} = \ldots - k^2 \mathcal F\left[\mathcal F^{-1}[i k \phi_t]^2\right] $$