Tutorial 03: Boundary conditions

The particle systems are defined in a rectangular box whose dimensions are specified by the attribute L.

The boundary conditions are specified by the attribute bc which can be one of the following.

• A list of size \(d\) containing for each dimension either 0 (periodic) or 1 (wall with reflecting boundary conditions).

• The string "open" : no boundary conditions.

• The string "periodic" : periodic boundary conditions.

• The string "spherical" : reflecting boundary conditions on the sphere of diameter \(L\) enclosed in the square domain \([0,L]^d\).

For instance, let us simulate the Vicsek model in an elongated rectangular domain \([0,L_x]\times[0,L_y]\) with periodic boundary conditions in the \(x\)-dimension and reflecting boundary conditions in the \(y\)-dimension.

First, some standard imports…

import time
import torch
from sisyphe.models import Vicsek
from sisyphe.display import display_kinetic_particles

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

The parameters of the model

N = 100000

R = .01
c = .1
nu = 3.
sigma = 1.

dt = .01

variant = {"name" : "max_kappa", "parameters" : {"kappa_max" : 10.}}

The spatial domain, the boundary conditions and the initial conditions…

Lx = 3.
Ly = 1./3.
L = [Lx, Ly]
bc = [0,1]

pos = torch.rand((N,2)).type(dtype)
pos[:,0] = L[0]*pos[:,0]
pos[:,1] = L[1]*pos[:,1]
vel = torch.randn(N,2).type(dtype)
vel = vel/torch.norm(vel,dim=1).reshape((N,1))

simu = Vicsek(
    pos = pos.detach().clone(),
    vel = vel.detach().clone(),
    v = c,
    sigma = sigma,
    nu = nu,
    interaction_radius = R,
    box_size = L,
    boundary_conditions=bc,
    dt = dt,
    variant = variant,
    block_sparse_reduction = True,
    number_of_cells = 100**2)

Finally run the simulation over 300 units of time….

# sphinx_gallery_thumbnail_number = 15

frames = [0., 2., 5., 10., 30., 42., 71., 100, 123, 141, 182, 203, 256, 272, 300]

s = time.time()
it, op = display_kinetic_particles(simu, frames, order=True, figsize=(8,3))
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%

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: 523.2031211853027 seconds
Average time per iteration: 0.01744010403951009 seconds

The simulation produces small clusters moving from left to right or from right to left. Each “step” in the order parameter corresponds to a collision between two clusters moving in opposite directions.