Skip to content
Snippets Groups Projects
Commit a8707c7c authored by Jonas Plewinski's avatar Jonas Plewinski
Browse files

resolved errors in demo notebook and integrated unittest for shift_matrix

parent 621a9954
No related branches found
No related tags found
1 merge request!44Added central moments functionality to lbmpy
Pipeline #28160 waiting for manual action
Stage: test
Hide whitespace changes
Inline Side-by-side
doc/notebooks/demo_central_moments.ipynb
+ 18
− 18
View file @ a8707c7c
%% Cell type:code id: tags:
``` python
# imports
from lbmpy.session import *
from lbmpy.moments import *
from lbmpy.forcemodels import Guo
from lbmpy.methods.abstractlbmethod import LbmCollisionRule
from collections import OrderedDict
from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments
from functions import get_shift_matrix, get_central_moments
from collections import OrderedDict
```
%% Cell type:markdown id: tags:
# Tutorial: Central Moments
In this tutorial, we will learn what central moments are, what advantages they have and how we can use them with lbmpy.
%% Cell type:markdown id: tags:
## 1) Theoretical Background
%% Cell type:markdown id: tags:
### What are Central Moments and why use them?
The most common collision operators are the SRT and MRT operators. Both of these operators are based on a second-order approximation of the Maxwell-Boltzmann distribution, where the relaxation is done in either single or multiple steps. However, both approaches have the disadvantage that at relaxation rates that are close to the limit, both methods become unstable. For these kinds of problems, other operators can be used, e.g. the entropic LB operator, which enforces an H-theorem on the lattice and is therefore unconditionally stable. In order to get this feature of unconditional stability, the relaxation time is modulated in dependence of the local entropy, which removes high frequencies from the flux fields. Nevertheless, the method that is presented in this notebook is another one, the Cascaded LBM (= CLBM), which uses central moments. It does not remove any frequencies, but instead tries to deal with the high frequencies with sufficient accuracy, such that these do not threaten the overall stability.
%% Cell type:markdown id: tags:
So CLBM has the goal to remove instabilities which can be traced back to an insufficient level of Galilean invariance (GI). GI describes the independence of physical processes from the speed of the reference system.
%% Cell type:markdown id: tags:
From tutorial 3 we know, that _lbmpy_ is using the following equation for moment based relaxations
$$K(f) = f - C^{-1}S\left(Cf - c^{(eq)}\right)$$
A few things are happening here:
* $Cf$ is transforming the distribution function $f$ into the moment space. The solution of this linear transformation is the moment-vector $c = Cf$.
* $c^{(eq)}$ defines the equilibrium in the moment space (alternatively $Cf^{(eq)}$, but it is more natural to define it directly in the moment space)
* In the moment space, the relaxation is executed by multiplying $S$ from left, to get $\Delta m = S \cdot \left(m - m^{(eq)}\right)$
* then an inverse transformation of the relaxed moments is done by $\Delta f = C^{-1} \Delta m$
* and in the end, the streaming is happening by $K(f) = f - \Delta f$
%% Cell type:markdown id: tags:
Now, in order to restore the Galilean invariance, we will use central moments. For that, we have to perform an additional transformation; we have to transform the raw moments into the central moment space. First, we define the continuous central moments:
%% Cell type:markdown id: tags:
$$\kappa_{\alpha \beta \gamma} = \int_{\mathbb{R}^3} (\xi_x - u_x)^\alpha \cdot (\xi_y - u_y)^\beta \cdot (\xi_z - u_z)^\gamma \cdot f(t, x, \xi) \,\, d\xi$$
%% Cell type:markdown id: tags:
where $u = \left( u_x, u_y, u_z \right)$ describes the macroscopic velocities. The discrete central moments are defined as
%% Cell type:markdown id: tags:
$$\kappa_{\alpha \beta \gamma} = \sum_{i = 0}^{n_v} [(c_i)_x - u_x]^\alpha \cdot [(c_i)_y - u_y]^\beta \cdot [(c_i)_z - u_z]^\gamma \cdot f_i \\
\qquad\qquad\quad = \sum_{i, j, k = -1}^{1} [(c_{ijk})_x - u_x]^\alpha \cdot [(c_{ijk})_y - u_y]^\beta \cdot [(c_{ijk})_z - u_z]^\gamma \cdot f_{ijk} $$
%% Cell type:markdown id: tags:
The main difference between the transformation to raw moments and the one to central moments is that the first is a linear transformation, while the second is a polynomial transformation. As a consequence, we cannot easily formulate this transformation as a simple matrix-vector multiplication; nevertheless, we can use the fast central moment transformation:
%% Cell type:markdown id: tags:
$$ \kappa_{ij\,\mid\,\gamma} = \sum_{k = -1}^{1} [(c_{ijk})_z - u_z]^\gamma \cdot f_{ijk} $$
$$ \kappa_{i\,\mid\,\beta \gamma} = \sum_{j = -1}^{1} [(c_{ijk})_y - u_y]^\beta \cdot \kappa_{ij\,\mid\,\gamma} $$
$$ \kappa_{\alpha \beta \gamma} = \sum_{i = -1}^{1} [(c_{ijk})_x - u_x]^\alpha \cdot \kappa_{i\,\mid\, \beta \gamma} $$
%% Cell type:markdown id: tags:
Now we can replace the transformation from the distribution function into the raw moment space with the fast central moment transformation into the central moment space, but because it is not a linear transformation, we are not able to write this as one equation as before, at least not directly. Therefore we have to do computations step by step and put the results together, such that we get the following algorithm:
%% Cell type:markdown id: tags:
* Transform distribution function $f_{ijk}$ with the fast central moment transformation into the central moments $\kappa_{\alpha \beta \gamma}$ $\left(f_{ijk} \Rightarrow \kappa_{\alpha\beta\gamma}\right)$
* Calculate the central equilibrium moments $\kappa_{i}^{(eq)}$ $\left(f_{ijk}^{(eq)} \Rightarrow \kappa_{\alpha\beta\gamma}^{(eq)}\right)$
* Relax the central moments $\kappa_i$ each with the relaxation frequency $\omega_i$ (diagonal elements of $S$) $\left(\kappa^{*} = S \left( \kappa_{\alpha\beta\gamma} - \kappa_{\alpha\beta\gamma}^{(eq)} \right)\right)$
* Transform the relaxed central moments with the inverse central moments transformation into the relaxed distribution function $f_{ijk}^{*}$ $\left(\kappa^{*} \Rightarrow f^{*}\right)$
* Do the update $\left(f(t+\Delta t, x+c_x \Delta t) = f(t, x) - f^{*}(t, x)\right)$ and streaming step
%% Cell type:markdown id: tags:
### Possible Solutions to Fix the Problem
%% Cell type:markdown id: tags:
However, there are two possibilities on how to solve the dilemma and write this again in one equation. One is to modify the above-described approach, and the other one uses features of lbmpy for a workaround:
1. Extract with the help of the fast central moment transformation a shift matrix, which is lower triangular, to transform the moment space into the central moment space.
2. Set up a method directly with the central moments. However, the problem with this approach is, that the transformation matrix which brings the distribution function into the moment space changes, and is not sparse anymore but full. As a consequence, the subexpression elimination is not able to reduce the operation count as much as in the case of the shift matrix (D2Q9, in D3Q27 case the subexpression elimination is not working anymore), which we will later see in an example.
Let us first look at the first possibility. We introduce a little modification to the approach described above. With the help of the fast central moment transformation, we can extract a shift matrix $N$, which transforms the raw moments into the central moment space via matrix multiplication. The shift matrix consists of the coefficients of the moments from the polynomials, which we get from the fast central moment transformation. Furthermore, with this, we can again write the collision equation in a form with matrix-matrix/vector multiplications and do not have to execute multiple things step after step.
%% Cell type:markdown id: tags:
$$K(f) = f - C^{-1}N^{-1}SN\left(Cf - c^{(eq)}\right)$$
%% Cell type:markdown id: tags:
## 2) Using Central Moments with *lbmpy* by Using the Shift Matrix $N$
%% Cell type:markdown id: tags:
First let us define variables which we later need, like the velocity $u$, the distribution function $f$, the relaxation rates $\omega$ and the viscosity $\rho$.
%% Cell type:code id: tags:
``` python
u = sp.Matrix(sp.symbols("u_:2")) # velocity
f = sp.Matrix(sp.symbols("f_:9")) # distribution function
ω = sp.Matrix(sp.symbols("omega_:9")) # relaxation rates
ρ = sp.symbols("rho") # viscosity
```
%% Cell type:markdown id: tags:
We will show the usage of the central moments in an easy setup, the simulation of a force-driven channel flow. First, we define a method with raw moments.
%% Cell type:code id: tags:
``` python
stencil = get_stencil("D2Q9", ordering='walberla')
force = [5e-5, 0]
F = Guo(force)
method = create_lb_method(stencil=stencil, method='mrt_raw', compressible=True,
equilibrium_order=4, force_model="Guo", force=force)
method
```
%% Output
<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f9b5d1abfd0>
<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f7cb10fb1c0>
%% Cell type:markdown id: tags:
For the sake of simplicity, we use an SRT approach and set the single relaxation rate to 1.9.
%% Cell type:code id: tags:
``` python
def modification_func(moment, eq, rr):
rr = 1.99
return moment, eq, rr
method = create_lb_method_from_existing(method, modification_func)
method
```
%% Output
<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f9b5d2683d0>
<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f7cb0f67fd0>
%% Cell type:markdown id: tags:
With the help of the method, we can extract specific vectors/matrices, which we will need to define our collision equation (transformation matrix to transform distribution function $f$ into the moment space, the discrete equilibrium values in the moment space). Further, we use the extracted moments and the defined stencil to get the second transformation matrix with `get_shift_matrix`, to get from the moment space to the central moment space.
%% Cell type:code id: tags:
``` python
raw_moments = method.moments
C = moment_matrix(raw_moments, stencil)
S = sp.diag(*method.relaxation_rates)
c_eq = sp.Matrix(method.moment_equilibrium_values)
N = get_shift_matrix(raw_moments, stencil)
N = shift_matrix(raw_moments, stencil)
N
```
%% Output
$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- u_{0} & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- u_{1} & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\u_{0}^{2} & - 2 u_{0} & 0 & 1 & 0 & 0 & 0 & 0 & 0\\u_{1}^{2} & 0 & - 2 u_{1} & 0 & 1 & 0 & 0 & 0 & 0\\u_{0} u_{1} & - u_{1} & - u_{0} & 0 & 0 & 1 & 0 & 0 & 0\\- u_{0}^{2} u_{1} & 2 u_{0} u_{1} & u_{0}^{2} & - u_{1} & 0 & - 2 u_{0} & 1 & 0 & 0\\- u_{0} u_{1}^{2} & u_{1}^{2} & 2 u_{0} u_{1} & 0 & - u_{0} & - 2 u_{1} & 0 & 1 & 0\\u_{0}^{2} u_{1}^{2} & - 2 u_{0} u_{1}^{2} & - 2 u_{0}^{2} u_{1} & u_{1}^{2} & u_{0}^{2} & 4 u_{0} u_{1} & - 2 u_{1} & - 2 u_{0} & 1\end{matrix}\right]$
⎡ 1 0 0 0 0 0 0 0 0⎤
⎢ ⎥
⎢ -u₀ 1 0 0 0 0 0 0 0⎥
⎢ ⎥
⎢ -u₁ 0 1 0 0 0 0 0 0⎥
⎢ ⎥
⎢ 2 ⎥
⎢ u₀ -2⋅u₀ 0 1 0 0 0 0 0⎥
⎢ ⎥
⎢ 2 ⎥
⎢ u₁ 0 -2⋅u₁ 0 1 0 0 0 0⎥
⎢ ⎥
⎢ u₀⋅u₁ -u₁ -u₀ 0 0 1 0 0 0⎥
⎢ ⎥
⎢ 2 2 ⎥
⎢-u₀ ⋅u₁ 2⋅u₀⋅u₁ u₀ -u₁ 0 -2⋅u₀ 1 0 0⎥
⎢ ⎥
⎢ 2 2 ⎥
⎢-u₀⋅u₁ u₁ 2⋅u₀⋅u₁ 0 -u₀ -2⋅u₁ 0 1 0⎥
⎢ ⎥
⎢ 2 2 2 2 2 2 ⎥
⎣u₀ ⋅u₁ -2⋅u₀⋅u₁ -2⋅u₀ ⋅u₁ u₁ u₀ 4⋅u₀⋅u₁ -2⋅u₁ -2⋅u₀ 1⎦
%% Cell type:markdown id: tags:
With the definition of all relevant variables, vectors and matrices, we can define the collision equation as described above.
%% Cell type:code id: tags:
``` python
collision_equation = f + C.inv() * N.inv() * S * N * (c_eq - C * f) + sp.Matrix(F(method))
```
%% Cell type:markdown id: tags:
In order to use this collision equation later on in a predefined scenario, we have to convert it into a collision rule, which we will do below, while introducing three subexpressions for the calculation of $\rho$, $u_0$ and $u_1$.
%% Cell type:code id: tags:
``` python
collision_rule = [Assignment(lhs, rhs) for lhs, rhs in zip(method.post_collision_pdf_symbols, collision_equation)]
@ps.kernel
def subexpressions():
ρ @= sum(f)
u[0] @= sum(f.T * sp.Matrix(np.array(stencil)[:, 0]))/ρ + force[0]/(2*ρ)
u[1] @= sum(f.T * sp.Matrix(np.array(stencil)[:, 1]))/ρ + force[1]/(2*ρ)
collision_rule = LbmCollisionRule(method, main_assignments=collision_rule, subexpressions=[subexpressions[0],
subexpressions[1],
subexpressions[2]])
```
%% Cell type:markdown id: tags:
After that, we can use the subexpression elimination from pystencils to reduce the number of operations needed to perform the collision.
%% Cell type:code id: tags:
``` python
generic_strategy = ps.simp.SimplificationStrategy()
generic_strategy.add(ps.simp.sympy_cse)
generic_strategy.create_simplification_report(collision_rule)
```
%% Output
<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f9b78f46610>
<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f7cb0f678b0>
%% Cell type:code id: tags:
``` python
collision_rule = ps.simp.sympy_cse(collision_rule)
```
%% Cell type:code id: tags:
``` python
collision_rule.operation_count
```
%% Output
{'adds': 393,
'muls': 246,
'divs': 1,
'sqrts': 0,
'fast_sqrts': 0,
'fast_inv_sqrts': 0,
'fast_div': 0}
%% Cell type:markdown id: tags:
Now we can pass our method and customized collision rule to the predefined scenario `create_channel` to simulate an easy channel flow.
%% Cell type:code id: tags:
``` python
ch_with_cm = create_channel((300,100), method=method, force=force[0], force_model="Guo", initial_velocity=(0.5, 0), collision_rule=collision_rule,
optimization={"target": "cpu", "openmp": 4})
ch_with_cm.run(10000)
plt.vector_field_magnitude(ch_with_cm.velocity[:, :])
```
%% Output
<matplotlib.image.AxesImage at 0x7f9b78828280>
<matplotlib.image.AxesImage at 0x7f7ccc65cac0>
%% Cell type:markdown id: tags:
Just for verification, we simulate the same setup, but with a standard SRT approach, without central moments.
%% Cell type:code id: tags:
``` python
ch_normal_srt = create_channel(domain_size=(300,100), force=force[0], initial_velocity=(0.5,0),
relaxation_rate=1.99, optimization={'target': 'cpu'})
ch_normal_srt.run(10000)
plt.vector_field_magnitude(ch_normal_srt.velocity[:, :])
```
%% Output
<matplotlib.image.AxesImage at 0x7f9b73f22d00>
<matplotlib.image.AxesImage at 0x7f7cc7f6d190>
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(13, 10))
ax = plt.gca()
fontsize = 22
vel_profile_ch_with_cm = ch_with_cm.velocity[0.5, :, 0]
vel_profile_ch_normal_srt = ch_normal_srt.velocity[0.5, :, 0]
plt.xlabel('y-Direction', fontsize=fontsize)
plt.ylabel('Velocity Magnitude', fontsize=fontsize)
plt.grid(color='black', linestyle='-', linewidth=0.1)
plt.rc('font', size=fontsize)
plt.plot(vel_profile_ch_normal_srt, linestyle='-.', color='r')
plt.plot(vel_profile_ch_with_cm, color='b', linewidth=0.7)
plt.legend(['Velocity Profile Standard SRT', 'Velocity Profile Central Moments'], fontsize=fontsize)
```
%% Output
<matplotlib.legend.Legend at 0x7f9b73ed9370>
<matplotlib.legend.Legend at 0x7f7cc7f18b80>
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(13, 10))
ax = plt.gca()
fontsize = 22
plt.xlabel('y-Direction', fontsize=fontsize)
plt.ylabel('Difference in Velocity Magnitude', fontsize=fontsize)
plt.grid(color='black', linestyle='-', linewidth=0.1)
plt.rc('font', size=fontsize)
plt.plot(vel_profile_ch_with_cm - vel_profile_ch_normal_srt, color='b')
```
%% Output
[<matplotlib.lines.Line2D at 0x7f9b73e689a0>]
[<matplotlib.lines.Line2D at 0x7f7cc7f45730>]
%% Cell type:markdown id: tags:
### Compare Operation Counts of Both Described Methods
%% Cell type:markdown id: tags:
First we set up a method directly with the central moments.
%% Cell type:code id: tags:
``` python
cm = get_central_moments(u = u)
cm = get_central_moments(stencil, u = u)
mtrr = OrderedDict(zip(cm, ω))
```
%% Cell type:code id: tags:
``` python
method_dir = create_with_discrete_maxwellian_eq_moments(stencil, mtrr, compressible=True, equilibrium_order=4)
def modification_func(moment, eq, rr):
rr = 1.99
return moment, eq.subs(u[0], 0).subs(u[1], 0), rr
method_dir = create_lb_method_from_existing(method_dir, modification_func)
method_dir
```
%% Output
<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f9b73948f70>
<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f7ccc63c850>
%% Cell type:markdown id: tags:
Now let us first look at the moment matrix, to show, that it is indeed not sparse anymore.
%% Cell type:code id: tags:
``` python
method_dir.moment_matrix
```
%% Output
$\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\- u_{0} & - u_{0} & - u_{0} & - u_{0} - 1 & 1 - u_{0} & - u_{0} - 1 & 1 - u_{0} & - u_{0} - 1 & 1 - u_{0}\\- u_{1} & 1 - u_{1} & - u_{1} - 1 & - u_{1} & - u_{1} & 1 - u_{1} & 1 - u_{1} & - u_{1} - 1 & - u_{1} - 1\\u_{0}^{2} & u_{0}^{2} & u_{0}^{2} & \left(- u_{0} - 1\right)^{2} & \left(1 - u_{0}\right)^{2} & \left(- u_{0} - 1\right)^{2} & \left(1 - u_{0}\right)^{2} & \left(- u_{0} - 1\right)^{2} & \left(1 - u_{0}\right)^{2}\\u_{1}^{2} & \left(1 - u_{1}\right)^{2} & \left(- u_{1} - 1\right)^{2} & u_{1}^{2} & u_{1}^{2} & \left(1 - u_{1}\right)^{2} & \left(1 - u_{1}\right)^{2} & \left(- u_{1} - 1\right)^{2} & \left(- u_{1} - 1\right)^{2}\\u_{0} u_{1} & - u_{0} \left(1 - u_{1}\right) & - u_{0} \left(- u_{1} - 1\right) & - u_{1} \left(- u_{0} - 1\right) & - u_{1} \left(1 - u_{0}\right) & \left(1 - u_{1}\right) \left(- u_{0} - 1\right) & \left(1 - u_{0}\right) \left(1 - u_{1}\right) & \left(- u_{0} - 1\right) \left(- u_{1} - 1\right) & \left(1 - u_{0}\right) \left(- u_{1} - 1\right)\\- u_{0}^{2} u_{1} & u_{0}^{2} \left(1 - u_{1}\right) & u_{0}^{2} \left(- u_{1} - 1\right) & - u_{1} \left(- u_{0} - 1\right)^{2} & - u_{1} \left(1 - u_{0}\right)^{2} & \left(1 - u_{1}\right) \left(- u_{0} - 1\right)^{2} & \left(1 - u_{0}\right)^{2} \left(1 - u_{1}\right) & \left(- u_{0} - 1\right)^{2} \left(- u_{1} - 1\right) & \left(1 - u_{0}\right)^{2} \left(- u_{1} - 1\right)\\- u_{0} u_{1}^{2} & - u_{0} \left(1 - u_{1}\right)^{2} & - u_{0} \left(- u_{1} - 1\right)^{2} & u_{1}^{2} \left(- u_{0} - 1\right) & u_{1}^{2} \left(1 - u_{0}\right) & \left(1 - u_{1}\right)^{2} \left(- u_{0} - 1\right) & \left(1 - u_{0}\right) \left(1 - u_{1}\right)^{2} & \left(- u_{0} - 1\right) \left(- u_{1} - 1\right)^{2} & \left(1 - u_{0}\right) \left(- u_{1} - 1\right)^{2}\\u_{0}^{2} u_{1}^{2} & u_{0}^{2} \left(1 - u_{1}\right)^{2} & u_{0}^{2} \left(- u_{1} - 1\right)^{2} & u_{1}^{2} \left(- u_{0} - 1\right)^{2} & u_{1}^{2} \left(1 - u_{0}\right)^{2} & \left(1 - u_{1}\right)^{2} \left(- u_{0} - 1\right)^{2} & \left(1 - u_{0}\right)^{2} \left(1 - u_{1}\right)^{2} & \left(- u_{0} - 1\right)^{2} \left(- u_{1} - 1\right)^{2} & \left(1 - u_{0}\right)^{2} \left(- u_{1} - 1\right)^{2}\end{matrix}\right]$
⎡ 1 1 1 1 1
⎢
⎢ -u₀ -u₀ -u₀ -u₀ - 1 1 - u₀
⎢
⎢ -u₁ 1 - u₁ -u₁ - 1 -u₁ -u₁
⎢
⎢ 2 2 2 2 2
⎢ u₀ u₀ u₀ (-u₀ - 1) (1 - u₀) (
⎢
⎢ 2 2 2 2 2
⎢ u₁ (1 - u₁) (-u₁ - 1) u₁ u₁ (
⎢
⎢ u₀⋅u₁ -u₀⋅(1 - u₁) -u₀⋅(-u₁ - 1) -u₁⋅(-u₀ - 1) -u₁⋅(1 - u₀) (1 -
⎢
⎢ 2 2 2 2 2
⎢-u₀ ⋅u₁ u₀ ⋅(1 - u₁) u₀ ⋅(-u₁ - 1) -u₁⋅(-u₀ - 1) -u₁⋅(1 - u₀) (1 - u
⎢
⎢ 2 2 2 2 2
⎢-u₀⋅u₁ -u₀⋅(1 - u₁) -u₀⋅(-u₁ - 1) u₁ ⋅(-u₀ - 1) u₁ ⋅(1 - u₀) (1 - u
⎢
⎢ 2 2 2 2 2 2 2 2 2 2
⎣u₀ ⋅u₁ u₀ ⋅(1 - u₁) u₀ ⋅(-u₁ - 1) u₁ ⋅(-u₀ - 1) u₁ ⋅(1 - u₀) (1 - u
1 1 1 1
-u₀ - 1 1 - u₀ -u₀ - 1 1 - u₀
1 - u₁ 1 - u₁ -u₁ - 1 -u₁ - 1
2 2 2 2
-u₀ - 1) (1 - u₀) (-u₀ - 1) (1 - u₀)
2 2 2 2
1 - u₁) (1 - u₁) (-u₁ - 1) (-u₁ - 1)
u₁)⋅(-u₀ - 1) (1 - u₀)⋅(1 - u₁) (-u₀ - 1)⋅(-u₁ - 1) (1 - u₀)⋅(-u₁ - 1
2 2 2 2
₁)⋅(-u₀ - 1) (1 - u₀) ⋅(1 - u₁) (-u₀ - 1) ⋅(-u₁ - 1) (1 - u₀) ⋅(-u₁ - 1
2 2 2
₁) ⋅(-u₀ - 1) (1 - u₀)⋅(1 - u₁) (-u₀ - 1)⋅(-u₁ - 1) (1 - u₀)⋅(-u₁ - 1)
2 2 2 2 2 2 2
₁) ⋅(-u₀ - 1) (1 - u₀) ⋅(1 - u₁) (-u₀ - 1) ⋅(-u₁ - 1) (1 - u₀) ⋅(-u₁ - 1
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
) ⎥
⎥
⎥
) ⎥
⎥
2 ⎥
⎥
⎥
2⎥
) ⎦
%% Cell type:markdown id: tags:
Next, we want to know how many operations the collision rule has and how much we can reduce them.
%% Cell type:code id: tags:
``` python
collision_rule_dir = method_dir.get_collision_rule()
```
%% Cell type:code id: tags:
``` python
generic_strategy_dir = ps.simp.SimplificationStrategy()
generic_strategy_dir.add(ps.simp.sympy_cse)
generic_strategy_dir.create_simplification_report(collision_rule_dir)
```
%% Output
<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f9b738c1100>
<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f7cc7993190>
%% Cell type:code id: tags:
``` python
collision_rule_dir = ps.simp.sympy_cse(collision_rule_dir)
```
%% Cell type:code id: tags:
``` python
collision_rule_dir.operation_count
```
%% Output
{'adds': 248,
'muls': 179,
'divs': 1,
'sqrts': 0,
'fast_sqrts': 0,
'fast_inv_sqrts': 0,
'fast_div': 0}
%% Cell type:code id: tags:
``` python
collision_rule.operation_count
```
%% Output
{'adds': 393,
'muls': 246,
'divs': 1,
'sqrts': 0,
'fast_sqrts': 0,
'fast_inv_sqrts': 0,
'fast_div': 0}
%% Cell type:markdown id: tags:
As we can see, not only is the initial operation count - before the subexpression elimination - much larger than the operation count with the shift matrix (116002 vs 2106 operations), but also the operation count after the common subexpression elimination is larger, and that is already the case for D2Q9, for D3Q27 the difference will get even more extensive.
%% Cell type:markdown id: tags:
## Literature
%% Cell type:markdown id: tags:
- Martin Geier, Andreas Greiner, and Jan G. Korvink. "Cascaded digital lattice Boltzmann automata for high Reynolds number flow". In: _Phys. Rev. E_ 73 (6. June 2006), p. 066705. DOI: 10.1103/PhysRevE.73.066705. URL: https://link.aps.prg/doi/10.1103/PhysRevE.73.066705.
- G. Gruszczyński et al. "A cascaded phase-field lattice Boltzmann model for the simulation of incompressible, immiscible fluids with high density contrast". In: _Computers Mathematics with Applications_ 79.4 (2020), pp. 1049-1071. ISSN:0898-1221. DOI: https://doi.org/10.1016/j.camwa.2019.08.018. URL: http://www.sciencedirect.com/science/article/pii/S0898122119304158
- Markus Muhr. "Kumulanten-basierte Lattice-Boltzmann-Methode". Seminar Ausarbeitung in _Transportprozesse: Mathematische Modellierung und Numerische Methoden_. Technische Universität München Fakultät für Mathematik.
......
lbmpy/moments.py
+ 3
− 4
View file @ a8707c7c
......@@ -379,14 +379,13 @@ def shift_matrix(moments, stencil):
directions = asarray(stencil)
u = sp.symbols("u_:{dim}".format(dim=dim))
f = sp.symbols("f_:{dim}".format(dim=nr_directions))
m = sp.symbols("m_:{dim}".format(dim=nr_directions))
u = sp.symbols(f"u_:{dim}")
f = sp.symbols(f"f_:{nr_directions}")
m = sp.symbols(f"m_:{nr_directions}")
Mf = moment_matrix(moments, stencil=stencil) * sp.Matrix(f)
shift_matrix_list = []
shift_equations = []
for nr_moment in range(len(moments)):
exponent_central_moment = moments[nr_moment].as_powers_dict()
......
lbmpy_tests/test_moments.py
+ 22
− 0
View file @ a8707c7c
from lbmpy.moments import *
from lbmpy.session import create_lb_method
from lbmpy.stencils import get_stencil
......@@ -100,3 +101,24 @@ def test_is_shear_moment():
assert is_shear_moment(x * y, 2)
assert is_shear_moment(x * y - 1, 2)
assert is_shear_moment(x * y - x, 2)
def test_shift_matrix():
stencils = (get_stencil("D2Q9", ordering='walberla'),
get_stencil("D3Q27", ordering='walberla'),
get_stencil("D3Q19", ordering='walberla'))
for stencil in stencils:
method = create_lb_method(stencil=stencil, method='mrt_raw', compressible=True, equilibrium_order=6)
moments = method.moments
sm = shift_matrix(moments, stencil)
m_eq = sp.Matrix(method.moment_equilibrium_values)
m_eq_rel = sp.simplify(sm * m_eq)
m_eq_sol = sp.Matrix(method.moment_equilibrium_values)
for i in range(0, len(stencil)):
m_eq_sol[i] = m_eq_sol[i].subs({sp.Symbol("u_0"):0, sp.Symbol("u_1"):0, sp.Symbol("u_2"):0})
assert m_eq_rel == m_eq_sol
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment