Skip to content
GitLab
Commit 488f8ef3 authored by Helen Schottenhamml's avatar Helen Schottenhamml
Browse files

Merge branch 'CentralMoments' into 'master' 

Central moments

See merge request pycodegen/lbmpy!71
parents 168d057a 86baa37b
Expand all Hide whitespace changes
Inline Side-by-side
doc/notebooks/demo_create_method_from_scratch.ipynb
View file @ 488f8ef3
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
```
%% Cell type:markdown id: tags:
# Demo: Create lbmpy Method from Scratch
<img src='../img/collision_space.svg' width="90%">
### Defining transformation to collision space
%% Cell type:code id: tags:
``` python
from lbmpy.moments import moment_matrix, moments_up_to_component_order, exponents_to_polynomial_representations
moment_exponents = list(moments_up_to_component_order(2, 2))
moment_exponents
```
%% Output
$\displaystyle \left[ \left( 0, \ 0\right), \ \left( 0, \ 1\right), \ \left( 0, \ 2\right), \ \left( 1, \ 0\right), \ \left( 1, \ 1\right), \ \left( 1, \ 2\right), \ \left( 2, \ 0\right), \ \left( 2, \ 1\right), \ \left( 2, \ 2\right)\right]$
[(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)]
%% Cell type:code id: tags:
``` python
moments = exponents_to_polynomial_representations(moment_exponents)
moments
```
%% Output
$\displaystyle \left( 1, \ y, \ y^{2}, \ x, \ x y, \ x y^{2}, \ x^{2}, \ x^{2} y, \ x^{2} y^{2}\right)$
⎛ 2 2 2 2 2 2⎞
⎝1, y, y , x, x⋅y, x⋅y , x , x ⋅y, x ⋅y ⎠
%% Cell type:code id: tags:
``` python
from lbmpy.stencils import get_stencil
d2q9 = get_stencil("D2Q9", ordering='walberla')
ps.stencil.plot(d2q9)
```
%% Output
%% Cell type:code id: tags:
``` python
M = moment_matrix(moments, stencil=d2q9)
M
```
%% Output
$\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1\\0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1\\0 & 0 & 0 & -1 & 1 & -1 & 1 & -1 & 1\\0 & 0 & 0 & 0 & 0 & -1 & 1 & 1 & -1\\0 & 0 & 0 & 0 & 0 & -1 & 1 & -1 & 1\\0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1\\0 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1\\0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\end{matrix}\right]$
⎡1 1 1 1 1 1 1 1 1 ⎤
⎢ ⎥
⎢0 1 -1 0 0 1 1 -1 -1⎥
⎢ ⎥
⎢0 1 1 0 0 1 1 1 1 ⎥
⎢ ⎥
⎢0 0 0 -1 1 -1 1 -1 1 ⎥
⎢ ⎥
⎢0 0 0 0 0 -1 1 1 -1⎥
⎢ ⎥
⎢0 0 0 0 0 -1 1 -1 1 ⎥
⎢ ⎥
⎢0 0 0 1 1 1 1 1 1 ⎥
⎢ ⎥
⎢0 0 0 0 0 1 1 -1 -1⎥
⎢ ⎥
⎣0 0 0 0 0 1 1 1 1 ⎦
%% Cell type:code id: tags:
``` python
from lbmpy.maxwellian_equilibrium import get_moments_of_continuous_maxwellian_equilibrium
from lbmpy.maxwellian_equilibrium import get_equilibrium_values_of_maxwell_boltzmann_function
eq_moments = get_moments_of_continuous_maxwellian_equilibrium(moments, order=2, dim=2,
c_s_sq=sp.Rational(1, 3))
eq_moments = get_equilibrium_values_of_maxwell_boltzmann_function(moments, order=2, dim=2,
c_s_sq=sp.Rational(1, 3),
space="moment")
omega = sp.symbols("omega")
relaxation_info = [(moment, eq_value, omega) for moment, eq_value in zip(moments, eq_moments)]
relaxation_info
```
%% Output
$\displaystyle \left[ \left( 1, \ \rho, \ \omega\right), \ \left( y, \ \rho u_{1}, \ \omega\right), \ \left( y^{2}, \ \rho u_{1}^{2} + \frac{\rho}{3}, \ \omega\right), \ \left( x, \ \rho u_{0}, \ \omega\right), \ \left( x y, \ \rho u_{0} u_{1}, \ \omega\right), \ \left( x y^{2}, \ \frac{\rho u_{0}}{3}, \ \omega\right), \ \left( x^{2}, \ \rho u_{0}^{2} + \frac{\rho}{3}, \ \omega\right), \ \left( x^{2} y, \ \frac{\rho u_{1}}{3}, \ \omega\right), \ \left( x^{2} y^{2}, \ \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}, \ \omega\right)\right]$
⎡
⎢ ⎛ 2 2 ρ ⎞
⎢(1, ρ, ω), (y, ρ⋅u₁, ω), ⎜y , ρ⋅u₁ + ─, ω⎟, (x, ρ⋅u₀, ω), (x⋅y, ρ⋅u₀⋅u₁, ω),
⎣ ⎝ 3 ⎠
⎛ 2 2
⎛ 2 ρ⋅u₀ ⎞ ⎛ 2 2 ρ ⎞ ⎛ 2 ρ⋅u₁ ⎞ ⎜ 2 2 ρ⋅u₀ ρ⋅u₁
⎜x⋅y , ────, ω⎟, ⎜x , ρ⋅u₀ + ─, ω⎟, ⎜x ⋅y, ────, ω⎟, ⎜x ⋅y , ───── + ───── +
⎝ 3 ⎠ ⎝ 3 ⎠ ⎝ 3 ⎠ ⎝ 3 3
⎞⎤
ρ ⎟⎥
─, ω⎟⎥
9 ⎠⎦
%% Cell type:code id: tags:
``` python
from lbmpy.methods.creationfunctions import create_generic_mrt
force_model = forcemodels.Guo(sp.symbols("F_:2"))
method = create_generic_mrt(d2q9, relaxation_info, compressible=False, force_model=force_model)
method
```
%% Output
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fc9d1ff5d60>
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f3aca401520>
%% Cell type:markdown id: tags:
### Example of a update equation without simplifications
%% Cell type:code id: tags:
``` python
collision_rule = method.get_collision_rule()
collision_rule
```
%% Output
AssignmentCollection: d_6, d_5, d_8, d_2, d_7, d_1, d_4, d_0, d_3 <- f(f_8, f_5, f_4, f_7, f_2, f_1, f_3, F_0, omega, F_1, f_0, f_6)
AssignmentCollection: d_1, d_2, d_7, d_4, d_3, d_0, d_8, d_5, d_6 <- f(omega, F_1, F_0, f_5, f_4, f_0, f_7, f_3, f_2, f_1, f_8, f_6)
%% Cell type:markdown id: tags:
### Generic simplification strategy - common subexpresssion elimination
%% 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 0x7fc9d1dd6610>
<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f3aca401190>
%% Cell type:markdown id: tags:
### A custom simplification strategy for moment-based methods
%% Cell type:code id: tags:
``` python
simplification_strategy = create_simplification_strategy(method)
simplification_strategy.create_simplification_report(collision_rule)
simplification_strategy.add(ps.simp.sympy_cse)
```
%% Cell type:markdown id: tags:
### Seeing the simplification in action
%% Cell type:code id: tags:
``` python
simplification_strategy.show_intermediate_results(collision_rule, symbols=[sp.Symbol("d_7")])
```
%% Output
<pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x7fc9d1da2880>
<pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x7f3aca39b160>
......
doc/sphinx/lbmpy.bib
View file @ 488f8ef3
......@@ -73,7 +73,6 @@ keywords = {lbm,multiphase,phasefield},
mendeley-tags = {lbm,multiphase,phasefield},
pages = {1--11},
title = {{Ternary free-energy lattice Boltzmann model with tunable surface tensions and contact angles}},
volume = {033305},
year = {2016}
}
......@@ -81,21 +80,27 @@ year = {2016}
author = {Geier, Martin and Sch{\"{o}}nherr, Martin and Pasquali, Andrea and Krafczyk, Manfred},
title = {{The cumulant lattice Boltzmann equation in three dimensions: Theory and validation}},
journal = {Computers \& Mathematics with Applications},
volume = {70},
number = {4},
pages = {507-547},
year = {2015},
doi = {10.1016/j.camwa.2015.05.001}
issn = {0898-1221},
doi = {10.1016/j.camwa.2015.05.001},
}
@Article{Coreixas2019,
author = {Christophe Coreixas and Bastien Chopard and Jonas Latt},
journal = {Physical Review E},
title = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
year = {2019},
month = {sep},
number = {3},
pages = {033305},
volume = {100},
doi = {10.1103/physreve.100.033305},
publisher = {American Physical Society ({APS})},
title = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
author = {Coreixas, Christophe and Chopard, Bastien and Latt, Jonas},
journal = {Phys. Rev. E},
volume = {100},
issue = {3},
pages = {033305},
numpages = {46},
year = {2019},
month = {Sep},
publisher = {American Physical Society},
doi = {10.1103/PhysRevE.100.033305},
url = {https://link.aps.org/doi/10.1103/PhysRevE.100.033305}
}
@PhdThesis{Geier2006,
......@@ -104,3 +109,14 @@ doi = {10.1016/j.camwa.2015.05.001}
title = {Ab inito derivation of the cascaded lattice Boltzmann automaton},
year = {2006},
}
@article{Fakhari2018,
title = {A phase-field lattice {Boltzmann} model for simulating multiphase flows in porous media: Application and comparison to experiments of {CO2} sequestration at pore scale},
journal = {Advances in Water Resources},
volume = {114},
pages = {119-134},
year = {2018},
issn = {0309-1708},
doi = {10.1016/j.advwatres.2018.02.005},
author = {Fakhari, A. and Li, Y. and Bolster, D. and Christensen, K. T.},
}
doc/sphinx/maxwellian_equilibrium.rst
View file @ 488f8ef3
......@@ -11,7 +11,7 @@ Maxwellian Equilibrium
.. autofunction:: lbmpy.maxwellian_equilibrium.continuous_maxwellian_equilibrium
.. autofunction:: lbmpy.maxwellian_equilibrium.get_moments_of_continuous_maxwellian_equilibrium
.. autofunction:: lbmpy.maxwellian_equilibrium.get_equilibrium_values_of_maxwell_boltzmann_function
.. autofunction:: lbmpy.maxwellian_equilibrium.get_moments_of_discrete_maxwellian_equilibrium
doc/sphinx/tutorials.rst
View file @ 488f8ef3
......@@ -17,6 +17,7 @@ You can open the notebooks directly to play around with the code examples.
/notebooks/07_tutorial_thermal_lbm.ipynb
/notebooks/08_tutorial_shanchen_twophase.ipynb
/notebooks/09_tutorial_shanchen_twocomponent.ipynb
/notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb
/notebooks/demo_stencils.ipynb
/notebooks/demo_create_method_from_scratch.ipynb
/notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb
......
lbmpy/continuous_distribution_measures.py
View file @ 488f8ef3
......@@ -6,11 +6,11 @@ import sympy as sp
from lbmpy.moments import polynomial_to_exponent_representation
from pystencils.cache import disk_cache, memorycache
from pystencils.sympyextensions import complete_the_squares_in_exp
from pystencils.sympyextensions import complete_the_squares_in_exp, scalar_product
@memorycache()
def moment_generating_function(generating_function, symbols, symbols_in_result):
def moment_generating_function(generating_function, symbols, symbols_in_result, velocity=None):
r"""
Computes the moment generating function of a probability distribution. It is defined as:
......@@ -21,6 +21,8 @@ def moment_generating_function(generating_function, symbols, symbols_in_result):
generating_function: sympy expression
symbols: a sequence of symbols forming the vector x
symbols_in_result: a sequence forming the vector t
velocity: if the generating function generates central moments, the velocity needs to be substracted. Thus the
velocity symbols need to be passed. All generating functions need to have the same parameters.
Returns:
transformation result F: an expression that depends now on symbols_in_result
......@@ -55,9 +57,27 @@ def moment_generating_function(generating_function, symbols, symbols_in_result):
return sp.simplify(result)
def cumulant_generating_function(func, symbols, symbols_in_result):
def central_moment_generating_function(func, symbols, symbols_in_result, velocity=sp.symbols("u_:3")):
r"""
Computes central moment generating func, which is defined as:
.. math ::
K( \vec{\Xi} ) = \exp ( - \vec{\Xi} \cdot \vec{u} ) M( \vec{\Xi}.
For parameter description see :func:`moment_generating_function`.
"""
Computes cumulant generating func, which is the logarithm of the moment generating func.
argument = - scalar_product(symbols_in_result, velocity)
return sp.exp(argument) * moment_generating_function(func, symbols, symbols_in_result)
def cumulant_generating_function(func, symbols, symbols_in_result, velocity=None):
r"""
Computes cumulant generating func, which is the logarithm of the moment generating func:
.. math ::
C(\vec{\Xi}) = \log M(\vec{\Xi})
For parameter description see :func:`moment_generating_function`.
"""
return sp.ln(moment_generating_function(func, symbols, symbols_in_result))
......@@ -93,16 +113,16 @@ def multi_differentiation(generating_function, index, symbols):
@memorycache(maxsize=512)
def __continuous_moment_or_cumulant(func, moment, symbols, generating_function):
def __continuous_moment_or_cumulant(func, moment, symbols, generating_function, velocity=sp.symbols("u_:3")):
if type(moment) is tuple and not symbols:
symbols = sp.symbols("xvar yvar zvar")
dim = len(moment) if type(moment) is tuple else len(symbols)
# not using sp.Dummy here - since it prohibits caching
t = tuple([sp.Symbol("tmpvar_%d" % i, ) for i in range(dim)])
t = sp.symbols(f"tmpvar_:{dim}")
symbols = symbols[:dim]
generating_function = generating_function(func, symbols, t)
generating_function = generating_function(func, symbols, t, velocity=velocity)
if type(moment) is tuple:
return multi_differentiation(generating_function, moment, t)
......@@ -128,6 +148,18 @@ def continuous_moment(func, moment, symbols=None):
return __continuous_moment_or_cumulant(func, moment, symbols, moment_generating_function)
def continuous_central_moment(func, moment, symbols=None, velocity=sp.symbols("u_:3")):
"""Computes central moment of given function.
Args:
func: function to compute moments of
moment: tuple or polynomial describing the moment
symbols: if moment is given as polynomial, pass the moment symbols, i.e. the dof of the polynomial
"""
return __continuous_moment_or_cumulant(func, moment, symbols, central_moment_generating_function,
velocity=velocity)
def continuous_cumulant(func, moment, symbols=None):
"""Computes cumulant of continuous function.
......
lbmpy/creationfunctions.py
View file @ 488f8ef3
......@@ -198,7 +198,8 @@ import lbmpy.forcemodels as forcemodels
import lbmpy.methods.centeredcumulant.force_model as cumulant_force_model
from lbmpy.fieldaccess import CollideOnlyInplaceAccessor, PdfFieldAccessor, PeriodicTwoFieldsAccessor
from lbmpy.fluctuatinglb import add_fluctuations_to_collision_rule
from lbmpy.methods import (create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc)
from lbmpy.methods import (create_mrt_orthogonal, create_mrt_raw, create_central_moment,
create_srt, create_trt, create_trt_kbc)
from lbmpy.methods.abstractlbmethod import RelaxationInfo
from lbmpy.methods.centeredcumulant import CenteredCumulantBasedLbMethod
from lbmpy.methods.momentbased.moment_transforms import PdfsToCentralMomentsByShiftMatrix
......@@ -417,7 +418,7 @@ def create_lb_method(**params):
'equilibrium_order': params['equilibrium_order'],
'force_model': force_model,
'maxwellian_moments': params['maxwellian_moments'],
'c_s_sq': params['c_s_sq'],
'c_s_sq': params['c_s_sq']
}
cumulant_params = {
......@@ -454,6 +455,8 @@ def create_lb_method(**params):
nested_moments = params['nested_moments'] if 'nested_moments' in params else None
method = create_mrt_orthogonal(stencil_entries, relaxation_rate_getter, weighted=weighted,
nested_moments=nested_moments, **common_params)
elif method_name.lower() == 'central_moment':
method = create_central_moment(stencil_entries, relaxation_rates, **common_params)
elif method_name.lower() == 'mrt_raw':
method = create_mrt_raw(stencil_entries, relaxation_rates, **common_params)
elif method_name.lower().startswith('trt-kbc-n'):
......@@ -529,7 +532,7 @@ def force_model_from_string(force_model_name, force_values):
'silva': forcemodels.Buick,
'edm': forcemodels.EDM,
'schiller': forcemodels.Schiller,
'cumulant': cumulant_force_model.CenteredCumulantForceModel
'cumulant': cumulant_force_model.CenteredCumulantForceModel,
}
if force_model_name.lower() not in force_model_dict:
raise ValueError("Unknown force model %s" % (force_model_name,))
......@@ -682,6 +685,6 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
stencil_entries = stencil_param
else:
stencil_entries = get_stencil(params_result['stencil'])
params_result['relaxation_rates'] = sp.symbols("omega_:%d" % len(stencil_entries))
params_result['relaxation_rates'] = sp.symbols(f"omega_:{len(stencil_entries)}")
return params_result, opt_params_result
lbmpy/cumulants.py
View file @ 488f8ef3
......@@ -124,7 +124,7 @@ def discrete_cumulant(func, cumulant, stencil):
assert len(stencil) == len(func)
dim = len(stencil[0])
wave_numbers = tuple([sp.Symbol("Xi_%d" % (i,)) for i in range(dim)])
wave_numbers = sp.symbols(f"Xi_:{dim}")
generating_function = __get_discrete_cumulant_generating_function(func, stencil, wave_numbers)
if type(cumulant) is tuple:
......
lbmpy/maxwellian_equilibrium.py
View file @ 488f8ef3
......@@ -10,6 +10,10 @@ import sympy as sp
from sympy import Rational as R
from pystencils.cache import disk_cache
from pystencils.sympyextensions import remove_higher_order_terms
from lbmpy.moments import MOMENT_SYMBOLS
from lbmpy.continuous_distribution_measures import continuous_moment, continuous_central_moment, continuous_cumulant
def get_weights(stencil, c_s_sq=sp.Rational(1, 3)):
......@@ -50,7 +54,7 @@ get_weights.weights = {
@disk_cache
def discrete_maxwellian_equilibrium(stencil, rho=sp.Symbol("rho"), u=tuple(sp.symbols("u_0 u_1 u_2")), order=2,
def discrete_maxwellian_equilibrium(stencil, rho=sp.Symbol("rho"), u=sp.symbols("u_:3"), order=2,
c_s_sq=sp.Symbol("c_s") ** 2, compressible=True):
"""
Returns the common discrete LBM equilibrium as a list of sympy expressions
......@@ -101,18 +105,19 @@ def discrete_maxwellian_equilibrium(stencil, rho=sp.Symbol("rho"), u=tuple(sp.sy
@disk_cache
def generate_equilibrium_by_matching_moments(stencil, moments, rho=sp.Symbol("rho"), u=tuple(sp.symbols("u_0 u_1 u_2")),
def generate_equilibrium_by_matching_moments(stencil, moments, rho=sp.Symbol("rho"), u=sp.symbols("u_:3"),
c_s_sq=sp.Symbol("c_s") ** 2, order=None):
"""
Computes discrete equilibrium, by setting the discrete moments to values taken from the continuous Maxwellian.
The number of moments has to match the number of directions in the stencil. For documentation of other parameters
see :func:`get_moments_of_continuous_maxwellian_equilibrium`
see :func:`get_equilibrium_values_of_maxwell_boltzmann_function`
"""
from lbmpy.moments import moment_matrix
dim = len(stencil[0])
Q = len(stencil)
assert len(moments) == Q, "Moment count(%d) does not match stencil size(%d)" % (len(moments), Q)
continuous_moments_vector = get_moments_of_continuous_maxwellian_equilibrium(moments, dim, rho, u, c_s_sq, order)
continuous_moments_vector = get_equilibrium_values_of_maxwell_boltzmann_function(moments, dim, rho, u, c_s_sq,
order, space="moment")
continuous_moments_vector = sp.Matrix(continuous_moments_vector)
M = moment_matrix(moments, stencil)
assert M.rank() == Q, "Rank of moment matrix (%d) does not match stencil size (%d)" % (M.rank(), Q)
......@@ -121,8 +126,8 @@ def generate_equilibrium_by_matching_moments(stencil, moments, rho=sp.Symbol("rh
@disk_cache
def continuous_maxwellian_equilibrium(dim=3, rho=sp.Symbol("rho"),
u=tuple(sp.symbols("u_0 u_1 u_2")),
v=tuple(sp.symbols("v_0 v_1 v_2")),
u=sp.symbols("u_:3"),
v=sp.symbols("v_:3"),
c_s_sq=sp.Symbol("c_s") ** 2):
"""
Returns sympy expression of Maxwell Boltzmann distribution
......@@ -141,15 +146,14 @@ def continuous_maxwellian_equilibrium(dim=3, rho=sp.Symbol("rho"),
return rho / (2 * sp.pi * c_s_sq) ** (sp.Rational(dim, 2)) * sp.exp(- vel_term / (2 * c_s_sq))
# -------------------------------- Equilibrium moments/cumulants ------------------------------------------------------
# -------------------------------- Equilibrium moments ----------------------------------------------------------------
@disk_cache
def get_moments_of_continuous_maxwellian_equilibrium(moments, dim, rho=sp.Symbol("rho"),
u=tuple(sp.symbols("u_0 u_1 u_2")),
c_s_sq=sp.Symbol("c_s") ** 2, order=None):
def get_equilibrium_values_of_maxwell_boltzmann_function(moments, dim, rho=sp.Symbol("rho"),
u=sp.symbols("u_:3"),
c_s_sq=sp.Symbol("c_s") ** 2, order=None,
space="moment"):
"""
Computes moments of the continuous Maxwell Boltzmann equilibrium distribution
Computes equilibrium values from the continuous Maxwell Boltzmann equilibrium.
Args:
moments: moments to compute, either in polynomial or exponent-tuple form
......@@ -159,19 +163,27 @@ def get_moments_of_continuous_maxwellian_equilibrium(moments, dim, rho=sp.Symbol
c_s_sq: symbol for speed of sound squared, defaults to symbol c_s**2
order: if this parameter is not None, terms that have a higher polynomial order in the macroscopic velocity
are removed
space: function space of the equilibrium values. Either moment, central moment or cumulant space are supported.
>>> get_moments_of_continuous_maxwellian_equilibrium( ( (0,0,0), (1,0,0), (0,1,0), (0,0,1), (2,0,0) ), dim=3 )
>>> get_equilibrium_values_of_maxwell_boltzmann_function( ( (0,0,0), (1,0,0), (0,1,0), (0,0,1), (2,0,0) ), dim=3 )
[rho, rho*u_0, rho*u_1, rho*u_2, rho*(c_s**2 + u_0**2)]
"""
from pystencils.sympyextensions import remove_higher_order_terms
from lbmpy.moments import MOMENT_SYMBOLS
from lbmpy.continuous_distribution_measures import continuous_moment
# trick to speed up sympy integration (otherwise it takes multiple minutes, or aborts):
# use a positive, real symbol to represent c_s_sq -> then replace this symbol afterwards with the real c_s_sq
c_s_sq_helper = sp.Symbol("csq_helper", positive=True, real=True)
mb = continuous_maxwellian_equilibrium(dim, rho, u, MOMENT_SYMBOLS[:dim], c_s_sq_helper)
result = [continuous_moment(mb, moment, MOMENT_SYMBOLS[:dim]).subs(c_s_sq_helper, c_s_sq) for moment in moments]
if space == "moment":
result = [continuous_moment(mb, moment, MOMENT_SYMBOLS[:dim]).subs(c_s_sq_helper, c_s_sq)
for moment in moments]
elif space == "central moment":
result = [continuous_central_moment(mb, moment, MOMENT_SYMBOLS[:dim], velocity=u).subs(c_s_sq_helper, c_s_sq)
for moment in moments]
elif space == "cumulant":
result = [continuous_cumulant(mb, moment, MOMENT_SYMBOLS[:dim]).subs(c_s_sq_helper, c_s_sq)
for moment in moments]
else:
raise ValueError("Only moment, central moment or cumulant space are supported")
if order is not None:
result = [remove_higher_order_terms(r, order=order, symbols=u) for r in result]
......@@ -180,7 +192,7 @@ def get_moments_of_continuous_maxwellian_equilibrium(moments, dim, rho=sp.Symbol
@disk_cache
def get_moments_of_discrete_maxwellian_equilibrium(stencil, moments,
rho=sp.Symbol("rho"), u=tuple(sp.symbols("u_0 u_1 u_2")),
rho=sp.Symbol("rho"), u=sp.symbols("u_:3"),
c_s_sq=sp.Symbol("c_s") ** 2, order=None, compressible=True):
"""Compute moments of discrete maxwellian equilibrium.
......@@ -235,32 +247,12 @@ def compressible_to_incompressible_moment_value(term, rho, u):
res += t
return res
# -------------------------------- Equilibrium moments -----------------------------------------------------------------
def get_cumulants_of_continuous_maxwellian_equilibrium(cumulants, dim, rho=sp.Symbol("rho"),
u=tuple(sp.symbols("u_0 u_1 u_2")), c_s_sq=sp.Symbol("c_s") ** 2,
order=None):
from lbmpy.moments import MOMENT_SYMBOLS
from lbmpy.continuous_distribution_measures import continuous_cumulant
from pystencils.sympyextensions import remove_higher_order_terms
# trick to speed up sympy integration (otherwise it takes multiple minutes, or aborts):
# use a positive, real symbol to represent c_s_sq -> then replace this symbol afterwards with the real c_s_sq
c_s_sq_helper = sp.Symbol("csq_helper", positive=True, real=True)
mb = continuous_maxwellian_equilibrium(dim, rho, u, MOMENT_SYMBOLS[:dim], c_s_sq_helper)
result = [continuous_cumulant(mb, cumulant, MOMENT_SYMBOLS[:dim]).subs(c_s_sq_helper, c_s_sq)
for cumulant in cumulants]
if order is not None:
result = [remove_higher_order_terms(r, order=order, symbols=u) for r in result]
return result
# -------------------------------- Equilibrium cumulants ---------------------------------------------------------------
@disk_cache
def get_cumulants_of_discrete_maxwellian_equilibrium(stencil, cumulants,
rho=sp.Symbol("rho"), u=tuple(sp.symbols("u_0 u_1 u_2")),
rho=sp.Symbol("rho"), u=sp.symbols("u_:3"),
c_s_sq=sp.Symbol("c_s") ** 2, order=None, compressible=True):
from lbmpy.cumulants import discrete_cumulant
if order is None:
......
lbmpy/methods/__init__.py
View file @ 488f8ef3
from lbmpy.methods.creationfunctions import (
create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc,
create_mrt_orthogonal, create_mrt_raw, create_central_moment, create_srt, create_trt, create_trt_kbc,
create_trt_with_magic_number, create_with_continuous_maxwellian_eq_moments,
create_with_discrete_maxwellian_eq_moments, mrt_orthogonal_modes_literature,
create_centered_cumulant_model, create_with_default_polynomial_cumulants,
......@@ -13,7 +13,7 @@ from .conservedquantitycomputation import DensityVelocityComputation
__all__ = ['RelaxationInfo', 'AbstractLbMethod',
'AbstractConservedQuantityComputation', 'DensityVelocityComputation',
'create_srt', 'create_trt', 'create_trt_with_magic_number', 'create_trt_kbc',
'create_mrt_orthogonal', 'create_mrt_raw',
'create_mrt_orthogonal', 'create_mrt_raw', 'create_central_moment',
'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments',
'mrt_orthogonal_modes_literature', 'create_centered_cumulant_model',
'create_with_default_polynomial_cumulants', 'create_with_polynomial_cumulants',
......
lbmpy/methods/abstractlbmethod.py
View file @ 488f8ef3
......@@ -32,12 +32,12 @@ class AbstractLbMethod(abc.ABC):
@property
def pre_collision_pdf_symbols(self):
"""Tuple of symbols representing the pdf values before collision"""
return sp.symbols("f_:%d" % (len(self.stencil),))
return sp.symbols(f"f_:{len(self.stencil)}")
@property
def post_collision_pdf_symbols(self):
"""Tuple of symbols representing the pdf values after collision"""
return sp.symbols("d_:%d" % (len(self.stencil),))
return sp.symbols(f"d_:{len(self.stencil)}")
# ------------------------- Abstract Methods & Properties ----------------------------------------------------------
......
lbmpy/methods/centeredcumulant/centered_cumulants.py
View file @ 488f8ef3
......@@ -6,10 +6,6 @@ from pystencils.stencil import have_same_entries
from lbmpy.moments import MOMENT_SYMBOLS, moment_sort_key, exponent_to_polynomial_representation
def statistical_quantity_symbol(name, exponents):
return sp.Symbol(f'{name}_{"".join(str(i) for i in exponents)}')
def exponent_tuple_sort_key(x):
return moment_sort_key(exponent_to_polynomial_representation(x))
......
lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
View file @ 488f8ef3
......@@ -14,12 +14,12 @@ from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantity
from lbmpy.moments import (
moments_up_to_order, get_order,
monomial_to_polynomial_transformation_matrix,
moment_sort_key, exponent_to_polynomial_representation, extract_monomials, MOMENT_SYMBOLS)
moment_sort_key, exponent_to_polynomial_representation, extract_monomials, MOMENT_SYMBOLS,
statistical_quantity_symbol)
# Local Imports
from lbmpy.methods.centeredcumulant.centered_cumulants import (
statistical_quantity_symbol, exponent_tuple_sort_key)
from lbmpy.methods.centeredcumulant.centered_cumulants import exponent_tuple_sort_key
from lbmpy.methods.centeredcumulant.cumulant_transform import (
PRE_COLLISION_CUMULANT, POST_COLLISION_CUMULANT,
......@@ -429,7 +429,7 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
if self._force_model is not None and \
not isinstance(self._force_model, CenteredCumulantForceModel) and include_force_terms:
force_model_terms = self._force_model(self)
force_term_symbols = sp.symbols("forceTerm_:%d" % (len(force_model_terms, )))
force_term_symbols = sp.symbols(f"forceTerm_:{len(force_model_terms)}")
force_subexpressions = [Assignment(sym, force_model_term)
for sym, force_model_term in zip(force_term_symbols, force_model_terms)]