# Solving a cable of cells on the GPU

In this demo we use the CUDA code generated by `gotranx` to solve a 1D cable of 500 cardiac
cells on an NVIDIA GPU, with one GPU thread per cell. A stimulus at one end of the cable starts
a wave of excitation that travels to the other end.

The cable is modelled with the monodomain equation, split into two steps per time step:

1. every cell's ODE (here the [ToR-ORd model](https://elifesciences.org/articles/48890)) is
   advanced with the [generalized Rush-Larsen scheme](../../docs/schemes/generalized_rush_larsen.md),
   in a kernel generated by `gotranx`;
2. the membrane potential diffuses to the neighbouring cells, in a small hand-written kernel.

```{note}
This page is not executed when the documentation is built, because that needs an NVIDIA GPU.
Its kernels are tested on the CPU, and the program is compiled with `nvcc` in the test suite.
```

## Generating the CUDA code

You can download the model in `.ode` format {download}`here <./ToRORd_dyn_chloride.ode>`. We
generate the kernels with

```shell
gotranx ode2cuda ToRORd_dyn_chloride.ode --scheme generalized_rush_larsen --namespace tor_ord -o tor_ord.cuh
```

which writes a header `tor_ord.cuh` with everything in the namespace `tor_ord`. Every function
of the model is a kernel that solves one cell per thread, and follows these conventions:

- Every kernel takes the number of cells, `n_cells`, as its last argument, and a thread whose
  cell is past the end returns at once. You can therefore launch more threads than cells.
- The states of all cells are stored **state-major**: state `i` of cell `c` is
  `states[i * n_cells + c]`. Neighbouring threads then read neighbouring memory, which is what
  a GPU is fastest at. It is the same layout as an array of shape `(num_states, n_cells)` in
  the Python backend.
- The parameters are shared by all cells.
- A scheme may update the states in place, i.e. you can pass the same array as `states` and
  `values`.
- `init_parameter_values` and the index lookups (`state_index`, `parameter_index`, ...) are
  ordinary host functions.
- The header contains kernels, so it can only be included from code compiled with `nvcc`
  (a `.cu` file), not from a `.cpp` file compiled with a host compiler.

For example, the scheme is

```c++
__global__ void generalized_rush_larsen(const double *states, const double t, const double dt,
                                        const double *__restrict__ parameters, double *values,
                                        const std::int64_t n_cells);
```

## The cable

Since the states are state-major, the potential of all cells is one contiguous row of the
states, `states + v * n_cells`, where `v = tor_ord::state_index("v")`. The kernel `diffuse` in
`cable.cuh` takes one explicit step of the diffusion of this row, with no flux through the
ends, and adds the stimulus to the first five cells during the first millisecond. `cable.cuh`
also holds the constants of the cable: a time step of 0.01 ms, 0.1 mm between cells and a
diffusion coefficient of 0.1 mm²/ms.

The program `main.cu` then

1. sets up the parameters on the host, switching off the model's own stimulus (it would
   stimulate every cell at once), and copies them to the GPU;
2. sets every cell's initial states with the kernel `tor_ord::init_state_values`;
3. for 500 ms, runs `tor_ord::generalized_rush_larsen` in place and then `cable::diffuse`;
4. every millisecond, copies the potential of all cells to the host and writes it as a row
   of `v.csv`.

## Running the demo

Generate `tor_ord.cuh` as above, then compile and run the program on a machine with an NVIDIA
GPU

```shell
nvcc -std=c++17 -O3 -arch=native -diag-suppress 177 main.cu -o cable
./cable v.csv
python3 plot.py v.csv
```

- `-arch=native` compiles for the GPU of the machine you compile on. If you compile on a
  machine without the GPU (e.g. the login node of a cluster), give the GPU's architecture
  instead, e.g. `-arch=sm_70` for a V100, `sm_80` for an A100 or `sm_90` for an H100. You can
  look it up with `nvidia-smi --query-gpu=name,compute_cap --format=csv` (7.0 means `sm_70`).
- `-diag-suppress 177` silences a warning about variables that are declared but not used,
  which the generated code has plenty of: every function assigns every parameter, whether it
  needs it or not.

`plot.py` draws the potential along the cable over time into `v.png`, next to `v.csv`:

```{figure} cable.png
:alt: The membrane potential along the cable over time

The membrane potential along the 50 mm cable over the first 500 ms, computed on a Tesla V100.
The wave of excitation starts at the stimulated end at 0 mm and reaches the far end after
about 90 ms. The cells repolarize about 360 ms after they were excited, so the wave of
repolarization travels along the cable in the same way.
```

The wave travels along the cable at about 0.55 m/s, which is within the range of the
conduction velocity in human ventricles.

## Troubleshooting

`main.cu` checks every CUDA call and stops with the error CUDA reports. The common ones are:

### `CUDA driver version is insufficient for CUDA runtime version`

The program cannot use the GPU's driver. Run `nvidia-smi`:

- If it says `No devices were found`, or there is no `nvidia-smi`, no GPU is visible to the
  program. On a cluster you are probably on a login node, or in a job without a GPU: run on a
  GPU node, e.g. with `srun --gres=gpu:1 ./cable v.csv` (depending on the cluster). A Docker
  container needs `--gpus all`.
- Otherwise, compare the `CUDA Version` in the top right of its output (the newest CUDA the
  driver supports) with `nvcc --version`. If `nvcc` is newer, e.g. CUDA 13.0 with a driver for
  CUDA 12.2, compile with an older CUDA (on a cluster often `module avail cuda`) or update the
  driver. Compiling with CUDA 12.x for a driver for an older CUDA 12 usually works, but CUDA 13
  needs a driver for CUDA 13.

### `no kernel image is available for execution on the device`

The program was compiled for another GPU architecture than the one it runs on. Compile with
`-arch=native` on the machine with the GPU, or with the architecture `nvidia-smi` reports (see
above).

```{note}
CUDA 13 no longer supports Volta GPUs (V100) and older. For those, compile with CUDA 12 and
`-arch=sm_70`.
```

## Source code

### `main.cu`
```{literalinclude} ./main.cu
:language: c++
```

### `cable.cuh`
```{literalinclude} ./cable.cuh
:language: c++
```

### `plot.py`
```{literalinclude} ./plot.py
:language: python
```