The storage

When you use pylbm, a generated code is performed using the descritpion of the scheme(s) (the velocities, the polynomials, the conserved moments, the equilibriums, …). There are several generators already implemented

To have best performance following the generator, you need a specific storage of the moments and distribution functions arrays. For example, it is preferable to have a storage like \([n_v, n_x, n_y, n_z]\) in NumPy \(n_v\) is the number of velocities and \(n_x\), \(n_y\) and \(n_z\) the grid size. It is due to the vectorized form of the algorithm. Whereas for Cython, it is preferable to have the storage \([n_x, n_y, n_z, n_v]\) using the pull algorithm.

So, we have implemented a storage class that always gives to the user the same access to the moments and disribution functions arrays but with a different storage in memory for the generator. This class is called Array.

It is really simple to create an array. You just need to give

2D example

Suppose that you want to create an array with a grid size \([5, 10]\) and \(9\) velocities with \(1\) cell in each direction for the fictitious domain.

[25]:
from pylbm.storage import Array
import numpy as np
a = Array(9, [5, 10], [1, 1])
[28]:
for i in range(a.nv):
    a[i] = i
[29]:
print(a[:])
[[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
  [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
  [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
  [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
  [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

 [[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
  [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
  [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
  [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
  [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]]

 [[ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
  [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
  [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
  [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
  [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]]

 [[ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]
  [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]
  [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]
  [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]
  [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]]

 [[ 4.  4.  4.  4.  4.  4.  4.  4.  4.  4.]
  [ 4.  4.  4.  4.  4.  4.  4.  4.  4.  4.]
  [ 4.  4.  4.  4.  4.  4.  4.  4.  4.  4.]
  [ 4.  4.  4.  4.  4.  4.  4.  4.  4.  4.]
  [ 4.  4.  4.  4.  4.  4.  4.  4.  4.  4.]]

 [[ 5.  5.  5.  5.  5.  5.  5.  5.  5.  5.]
  [ 5.  5.  5.  5.  5.  5.  5.  5.  5.  5.]
  [ 5.  5.  5.  5.  5.  5.  5.  5.  5.  5.]
  [ 5.  5.  5.  5.  5.  5.  5.  5.  5.  5.]
  [ 5.  5.  5.  5.  5.  5.  5.  5.  5.  5.]]

 [[ 6.  6.  6.  6.  6.  6.  6.  6.  6.  6.]
  [ 6.  6.  6.  6.  6.  6.  6.  6.  6.  6.]
  [ 6.  6.  6.  6.  6.  6.  6.  6.  6.  6.]
  [ 6.  6.  6.  6.  6.  6.  6.  6.  6.  6.]
  [ 6.  6.  6.  6.  6.  6.  6.  6.  6.  6.]]

 [[ 7.  7.  7.  7.  7.  7.  7.  7.  7.  7.]
  [ 7.  7.  7.  7.  7.  7.  7.  7.  7.  7.]
  [ 7.  7.  7.  7.  7.  7.  7.  7.  7.  7.]
  [ 7.  7.  7.  7.  7.  7.  7.  7.  7.  7.]
  [ 7.  7.  7.  7.  7.  7.  7.  7.  7.  7.]]

 [[ 8.  8.  8.  8.  8.  8.  8.  8.  8.  8.]
  [ 8.  8.  8.  8.  8.  8.  8.  8.  8.  8.]
  [ 8.  8.  8.  8.  8.  8.  8.  8.  8.  8.]
  [ 8.  8.  8.  8.  8.  8.  8.  8.  8.  8.]
  [ 8.  8.  8.  8.  8.  8.  8.  8.  8.  8.]]]
[30]:
b = Array(9, [5, 10], [1, 1], sorder=[2, 1, 0])
for i in range(b.nv):
    b[i] = i
[31]:
print(b[:])
[[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
  [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
  [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
  [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
  [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

 [[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
  [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
  [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
  [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
  [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]]

 [[ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
  [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
  [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
  [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
  [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]]

 [[ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]
  [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]
  [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]
  [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]
  [ 3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]]

 [[ 4.  4.  4.  4.  4.  4.  4.  4.  4.  4.]
  [ 4.  4.  4.  4.  4.  4.  4.  4.  4.  4.]
  [ 4.  4.  4.  4.  4.  4.  4.  4.  4.  4.]
  [ 4.  4.  4.  4.  4.  4.  4.  4.  4.  4.]
  [ 4.  4.  4.  4.  4.  4.  4.  4.  4.  4.]]

 [[ 5.  5.  5.  5.  5.  5.  5.  5.  5.  5.]
  [ 5.  5.  5.  5.  5.  5.  5.  5.  5.  5.]
  [ 5.  5.  5.  5.  5.  5.  5.  5.  5.  5.]
  [ 5.  5.  5.  5.  5.  5.  5.  5.  5.  5.]
  [ 5.  5.  5.  5.  5.  5.  5.  5.  5.  5.]]

 [[ 6.  6.  6.  6.  6.  6.  6.  6.  6.  6.]
  [ 6.  6.  6.  6.  6.  6.  6.  6.  6.  6.]
  [ 6.  6.  6.  6.  6.  6.  6.  6.  6.  6.]
  [ 6.  6.  6.  6.  6.  6.  6.  6.  6.  6.]
  [ 6.  6.  6.  6.  6.  6.  6.  6.  6.  6.]]

 [[ 7.  7.  7.  7.  7.  7.  7.  7.  7.  7.]
  [ 7.  7.  7.  7.  7.  7.  7.  7.  7.  7.]
  [ 7.  7.  7.  7.  7.  7.  7.  7.  7.  7.]
  [ 7.  7.  7.  7.  7.  7.  7.  7.  7.  7.]
  [ 7.  7.  7.  7.  7.  7.  7.  7.  7.  7.]]

 [[ 8.  8.  8.  8.  8.  8.  8.  8.  8.  8.]
  [ 8.  8.  8.  8.  8.  8.  8.  8.  8.  8.]
  [ 8.  8.  8.  8.  8.  8.  8.  8.  8.  8.]
  [ 8.  8.  8.  8.  8.  8.  8.  8.  8.  8.]
  [ 8.  8.  8.  8.  8.  8.  8.  8.  8.  8.]]]

You can see that the access of the data is the same for \(a\) et \(b\) whereas the sorder is not the same.

If we look at the array attribute which is the real storage of our data

[32]:
a.array.shape
[32]:
(9, 5, 10)
[33]:
b.array.shape
[33]:
(10, 5, 9)

you can see that it is not the same and it is exactly what we want. To do that, we use the swapaxes of numpy and we use this representation to have an access to our data.

Access to the data with the conserved moments

When you discribe your scheme, you define the conserved moments. It is usefull to have a direct acces to these moments by giving their name and not their indices in the array. So, it is possible to specify where are the conserved moments in the array.

Let define conserved moments using sympy symbol.

[35]:
import sympy
rho, u, v = sympy.symbols("rho, u, v")

We indicate to pylbm where are located these conserved moments in our array by giving a list of two elements: the first one is the scheme number and the second one the index in this scheme.

[45]:
a.set_conserved_moments({rho: [0, 0], u: [0, 2], v: [0, 1]})
[46]:
a[rho]
[46]:
array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])
[47]:
a[u]
[47]:
array([[ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.],
       [ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.],
       [ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.],
       [ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.],
       [ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.]])
[48]:
a[v]
[48]:
array([[ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]])
[ ]: