Note

Click here to download the full example code

Volume exclusion

This model is introduced in

Motsch, D. Peurichard, From short-range repulsion to Hele-Shaw problem in a model of tumor growth, J. Math. Biology, Vol. 76, No. 1, 2017.

First of all, some standard imports.

import os
import sys
import time
import math
import torch
import numpy as np
from matplotlib import pyplot as plt
import sisyphe.models as models
from sisyphe.display import scatter_particles


use_cuda = torch.cuda.is_available()
dtype = torch.cuda.FloatTensor if use_cuda else torch.FloatTensor

Repulsion force

Each particle is a disk with a (fixed) random radius. The particles repel each other when they overlap. The force exerted by a particle located at \(x_j\) with radius \(R_j\) on a particle located at \(x_i\) with radius \(R_i\) is

\[F = -\frac{\alpha}{R_i} \nabla_{x_i} U\left(\frac{|x_i - x_j|^2}{(R_i + R_j)^2}\right),\]

where the potential is

\[U(s) = -\log(s) + s - 1\,\,\text{for}\,\, s<1 \,\,\text{and}\,\, U(s) = 0\,\, \text{for}\,\, s>1.\]

Initially, the particles are clustered in a small region with a strong overlapping.

N = 10000
rmin = .1
rmax = 1.
R = (rmax-rmin)*torch.rand(N).type(dtype)+rmin
L = 100.
D0 = 20.
pos = (D0*torch.rand((N,2)).type(dtype)-D0/2)+torch.tensor([L/2,L/2]).type(dtype)

dt = .1

simu = models.VolumeExclusion(pos=pos,
                 interaction_radius=R,
                 box_size=L,
                 alpha=2.5,
                 division_rate=0.,
                 death_rate=0.,
                 dt=dt)

Run the simulation over 200 units of time using an adaptive time-step which ensures that the energy E of the system decreases.

# sphinx_gallery_thumbnail_number = 13

frames = [0,1,2,3,4,5,10,30,50,75,100,150,200]

s = time.time()
scatter_particles(simu,frames)
e = time.time()

Out:

Progress:0%
Progress:1%
Progress:2%
Progress:3%
Progress:4%
Progress:5%
Progress:6%
Progress:7%
Progress:8%
Progress:9%
Progress:10%
Progress:11%
Progress:12%
Progress:13%
Progress:14%
Progress:15%
Progress:16%
Progress:17%
Progress:18%
Progress:19%
Progress:20%
Progress:21%
Progress:22%
Progress:23%
Progress:24%
Progress:25%
Progress:26%
Progress:27%
Progress:28%
Progress:29%
Progress:30%
Progress:31%
Progress:32%
Progress:33%
Progress:34%
Progress:35%
Progress:36%
Progress:37%
Progress:38%
Progress:39%
Progress:40%
Progress:41%
Progress:42%
Progress:43%
Progress:44%
Progress:45%
Progress:46%
Progress:47%
Progress:48%
Progress:49%
Progress:50%
Progress:51%
Progress:52%
Progress:53%
Progress:54%
Progress:55%
Progress:56%
Progress:57%
Progress:58%
Progress:59%
Progress:60%
Progress:61%
Progress:62%
Progress:63%
Progress:64%
Progress:65%
Progress:66%
Progress:67%
Progress:68%
Progress:69%
Progress:70%
Progress:71%
Progress:72%
Progress:73%
Progress:74%
Progress:75%
Progress:76%
Progress:77%
Progress:78%
Progress:79%
Progress:80%
Progress:81%
Progress:82%
Progress:83%
Progress:84%
Progress:85%
Progress:86%
Progress:87%
Progress:88%
Progress:89%
Progress:90%
Progress:91%
Progress:92%
Progress:93%
Progress:94%
Progress:95%
Progress:96%
Progress:97%
Progress:98%
Progress:99%
Progress:100%

Print the total simulation time and the average time per iteration.

print('Total time: '+str(e-s)+' seconds')
print('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')

Out:

Total time: 326.2268440723419 seconds
Average time per iteration: 0.01703801347847401 seconds

Note this funny behaviour: the particles are clustered by size!

Repulsion force, random births and random deaths

Same system but this time, particles die at a constant rate and give birth to new particles at the same rate. A new particle is added next to its parent and has the same radius.

N = 10000
rmin = .1
rmax = 1.
R = (rmax-rmin)*torch.rand(N).type(dtype)+rmin
L = 100.
D0 = 20.
pos = (D0*torch.rand((N,2)).type(dtype)-D0/2)+torch.tensor([L/2,L/2]).type(dtype)

dt = .1

simu = models.VolumeExclusion(pos=pos,
                 interaction_radius=R,
                 box_size=L,
                 alpha=2.5,
                 division_rate=.3,
                 death_rate=.3,
                 dt=dt,
                 Nmax = 20000)

Run the simulation over 200 units of time using an adaptive time-step which ensures that the energy E of the system decreases.

frames = [0,1,2,3,4,5,10,30,50,75,100,150,200]

s = time.time()
scatter_particles(simu,frames)
e = time.time()

Out:

Progress:0%
Progress:1%
Progress:2%
Progress:3%
Progress:4%
Progress:5%
Progress:6%
Progress:7%
Progress:8%
Progress:9%
Progress:10%
Progress:11%
Progress:12%
Progress:13%
Progress:14%
Progress:15%
Progress:16%
Progress:17%
Progress:18%
Progress:19%
Progress:20%
Progress:21%
Progress:22%
Progress:23%
Progress:24%
Progress:25%
Progress:26%
Progress:27%
Progress:28%
Progress:29%
Progress:30%
Progress:31%
Progress:32%
Progress:33%
Progress:34%
Progress:35%
Progress:36%
Progress:37%
Progress:38%
Progress:39%
Progress:40%
Progress:41%
Progress:42%
Progress:43%
Progress:44%
Progress:45%
Progress:46%
Progress:47%
Progress:48%
Progress:49%
Progress:50%
Progress:51%
Progress:52%
Progress:53%
Progress:54%
Progress:55%
Progress:56%
Progress:57%
Progress:58%
Progress:59%
Progress:60%
Progress:61%
Progress:62%
Progress:63%
Progress:64%
Progress:65%
Progress:66%
Progress:67%
Progress:68%
Progress:69%
Progress:70%
Progress:71%
Progress:72%
Progress:73%
Progress:74%
Progress:75%
Progress:76%
Progress:77%
Progress:78%
Progress:79%
Progress:80%
Progress:81%
Progress:82%
Progress:83%
Progress:84%
Progress:85%
Progress:86%
Progress:87%
Progress:88%
Progress:89%
Progress:90%
Progress:91%
Progress:92%
Progress:93%
Progress:94%
Progress:95%
Progress:96%
Progress:97%
Progress:98%
Progress:99%
Progress:100%

Print the total simulation time and the average time per iteration.

print('Total time: '+str(e-s)+' seconds')
print('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')

Out:

Total time: 200.29545760154724 seconds
Average time per iteration: 0.007975450250917705 seconds

Total running time of the script: ( 8 minutes 53.860 seconds)

Gallery generated by Sphinx-Gallery