Skip to content
Snippets Groups Projects
Commit d77a9064 authored by itischler's avatar itischler
Browse files

Merge branch 'pressure_tensor_test' into 'pressure_tensor' 

Pressure tensor test

See merge request !1
parents 34272d9e 24239a5f
No related branches found
No related tags found
1 merge request!1Pressure tensor test
Pipeline #33245 failed
Stage: pretest
Stage: test
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1296 additions and 442 deletions
.gitlab-ci.yml
+ 31
− 29
View file @ d77a9064
stages:
- pretest
- test
- deploy
# -------------------------- Tests ------------------------------------------------------------------------------------
# -------------------------- Pre Tests --------------------------------------------------------------------------------
# Normal test - runs on every commit all but "long run" tests
tests-and-coverage:
stage: test
stage: pretest
except:
variables:
- $ENABLE_NIGHTLY_BUILDS
......@@ -36,6 +37,33 @@ tests-and-coverage:
cobertura: coverage.xml
junit: report.xml
minimal-conda:
stage: pretest
except:
variables:
- $ENABLE_NIGHTLY_BUILDS
image: i10git.cs.fau.de:5005/pycodegen/pycodegen/minimal_conda
script:
- pip install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
- python setup.py quicktest
tags:
- docker
# Linter for code formatting
flake8-lint:
stage: pretest
except:
variables:
- $ENABLE_NIGHTLY_BUILDS
image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
script:
- flake8 lbmpy
tags:
- docker
- cuda11
# -------------------------- Tests -------------------------------------------------------------------------------------
# pipeline with latest python version
latest-python:
stage: test
......@@ -132,18 +160,6 @@ ubuntu:
reports:
junit: report.xml
minimal-conda:
stage: test
except:
variables:
- $ENABLE_NIGHTLY_BUILDS
image: i10git.cs.fau.de:5005/pycodegen/pycodegen/minimal_conda
script:
- pip install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
- python setup.py quicktest
tags:
- docker
pycodegen-integration:
image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
stage: test
......@@ -179,21 +195,7 @@ pycodegen-integration:
- cuda11
- AVX
# -------------------- Linter & Documentation --------------------------------------------------------------------------
flake8-lint:
stage: test
except:
variables:
- $ENABLE_NIGHTLY_BUILDS
image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
script:
- flake8 lbmpy
tags:
- docker
- cuda11
# -------------------- Documentation and deploy ------------------------------------------------------------------------
build-documentation:
stage: test
......
AUTHORS.txt 0 → 100644
+ 13
− 0
View file @ d77a9064
Contributors:
-------------
- Martin Bauer <martin.bauer@fau.de>
- Markus Holzer <markus.holzer@fau.de>
- Michael Kuron <mkuron@icp.uni-stuttgart.de>
- Stephan Seitz <stephan.seitz@fau.de>
- Frederik Hennig <frederik.hennig@fau.de>
- Helen Schottenhamml <helen.schottenhamml@fau.de>
- Rudolf Weeber <weeber@icp.uni-stuttgart.de>
- Christian Godenschwager <christian.godenschwager@fau.de>
- Jan Hönig <jan.hoenig@fau.de>
CONTRIBUTING.md 0 → 100644
+ 3
− 0
View file @ d77a9064
# Contributing
lbmpy is built on the open-source python framework [pystencils](https://pypi.org/project/pystencils/). Please consider the [contribution guideline](https://i10git.cs.fau.de/pycodegen/pystencils/-/blob/master/CONTRIBUTING.md) of pystencils for contributing to lbmpy.
\ No newline at end of file
MANIFEST.in
+ 2
− 1
View file @ d77a9064
include README.md
include COPYING.txt
include RELEASE-VERSION
include AUTHORS.txt
include CONTRIBUTING.md
global-include *.pyx
include versioneer.py
include lbmpy/_version.py
README.md
+ 15
− 0
View file @ d77a9064
......@@ -55,3 +55,18 @@ Documentation
Read the docs [here](http://pycodegen.pages.i10git.cs.fau.de/lbmpy) and
check out the Jupyter notebooks in `doc/notebooks`.
Authors
-------
Many thanks go to the [contributors](AUTHORS.txt) of lbmpy.
### Please cite us
If you use lbmpy in a publication, please cite the following articles:
Overview:
- M. Bauer et al, lbmpy: Automatic code generation for efficient parallel lattice Boltzmann methods. Journal of Computational Science, 2021. https://doi.org/10.1016/j.jocs.2020.101269
Multiphase:
- M. Holzer et al, Highly efficient lattice Boltzmann multiphase simulations of immiscible fluids at high-density ratios on CPUs and GPUs through code generation. The International Journal of High Performance Computing Applications, 2021. https://doi.org/10.1177/10943420211016525
doc/notebooks/03_tutorial_lbm_formulation.ipynb
+ 166
− 21
View file @ d77a9064
This diff is collapsed.
doc/notebooks/04_tutorial_cumulant_LBM.ipynb
+ 48
− 272
View file @ d77a9064
This diff is collapsed.
doc/notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb 0 → 100644
+ 832
− 0
View file @ d77a9064
This diff is collapsed.
doc/notebooks/demo_create_method_from_scratch.ipynb
+ 9
− 8
View file @ d77a9064
%% 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
+ 28
− 12
View file @ d77a9064
......@@ -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
+ 1
− 1
View file @ d77a9064
......@@ -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
+ 1
− 0
View file @ d77a9064
......@@ -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/boundaries/__init__.py
+ 4
− 2
View file @ d77a9064
from lbmpy.boundaries.boundaryconditions import (
UBB, FixedDensity, SimpleExtrapolationOutflow, ExtrapolationOutflow, NeumannByCopy, NoSlip, StreamInConstant)
UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow,
ExtrapolationOutflow, NeumannByCopy, NoSlip, StreamInConstant)
from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
__all__ = ['NoSlip', 'UBB', 'SimpleExtrapolationOutflow', 'ExtrapolationOutflow', 'FixedDensity', 'NeumannByCopy',
__all__ = ['NoSlip', 'UBB', 'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
'FixedDensity', 'DiffusionDirichlet', 'NeumannByCopy',
'LatticeBoltzmannBoundaryHandling', 'StreamInConstant']
lbmpy/boundaries/boundaryconditions.py
+ 55
− 35
View file @ d77a9064
......@@ -79,6 +79,14 @@ class LbBoundary:
def name(self, new_value):
self._name = new_value
def __hash__(self):
return hash(self.name)
def __eq__(self, other):
if not isinstance(other, type(self)):
return False
return self.__dict__ == other.__dict__
# end class Boundary
......@@ -99,15 +107,6 @@ class NoSlip(LbBoundary):
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
return Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol))
def __hash__(self):
return hash(self.name)
def __eq__(self, other):
if not isinstance(other, NoSlip):
return False
return self.name == other.name
# end class NoSlip
......@@ -126,7 +125,6 @@ class UBB(LbBoundary):
"""
def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None, data_type='double'):
super(UBB, self).__init__(name)
self._velocity = velocity
self._adaptVelocityToForce = adapt_velocity_to_force
if callable(self._velocity) and not dim:
......@@ -136,12 +134,14 @@ class UBB(LbBoundary):
self.dim = dim
self.data_type = data_type
super(UBB, self).__init__(name)
@property
def additional_data(self):
""" In case of the UBB boundary additional data is a velocity vector. This vector is added to each cell to
realize velocity profiles for the inlet."""
if self.velocity_is_callable:
return [('vel_%d' % (i,), create_type(self.data_type)) for i in range(self.dim)]
return [(f'vel_{i}', create_type(self.data_type)) for i in range(self.dim)]
else:
return []
......@@ -212,7 +212,6 @@ class UBB(LbBoundary):
return [Assignment(f_in(inv_dir[direction]),
f_out(direction) - vel_term)]
# end class UBB
......@@ -259,6 +258,7 @@ class SimpleExtrapolationOutflow(LbBoundary):
tangential_offset = tuple(offset - normal for offset, normal in zip(neighbor_offset, self.normal_direction))
return Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](inv_dir[dir_symbol]))
# end class SimpleExtrapolationOutflow
......@@ -294,11 +294,10 @@ class ExtrapolationOutflow(LbBoundary):
streaming_pattern='pull', zeroth_timestep=Timestep.BOTH,
initial_density=None, initial_velocity=None, data_type='double'):
self.data_type = data_type
self.lb_method = lb_method
self.stencil = lb_method.stencil
self.dim = len(self.stencil[0])
if isinstance(normal_direction, str):
normal_direction = direction_string_to_offset(normal_direction, dim=self.dim)
......@@ -316,6 +315,8 @@ class ExtrapolationOutflow(LbBoundary):
self.initial_velocity = initial_velocity
self.equilibrium_calculation = None
self.data_type = data_type
if initial_density and initial_velocity:
equilibrium = lb_method.get_equilibrium(conserved_quantity_equations=AssignmentCollection([]))
rho = lb_method.zeroth_order_equilibrium_moment_symbol
......@@ -405,7 +406,6 @@ class ExtrapolationOutflow(LbBoundary):
return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
# end class ExtrapolationOutflow
......@@ -420,8 +420,9 @@ class FixedDensity(LbBoundary):
def __init__(self, density, name=None):
if name is None:
name = "Fixed Density " + str(density)
self.density = density
super(FixedDensity, self).__init__(name)
self._density = density
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
def remove_asymmetric_part_of_main_assignments(assignment_collection, degrees_of_freedom):
......@@ -441,7 +442,7 @@ class FixedDensity(LbBoundary):
density_symbol = cqc.defined_symbols()['density']
density = self._density
density = self.density
equilibrium_input = cqc.equilibrium_input_equations_from_init_values(density=density)
equilibrium_input = equilibrium_input.new_without_subexpressions()
density_eq = equilibrium_input.main_assignments[0]
......@@ -458,10 +459,43 @@ class FixedDensity(LbBoundary):
return subexpressions + [Assignment(f_in(inv_dir[dir_symbol]),
2 * eq_component - f_out(dir_symbol))]
# end class FixedDensity
class DiffusionDirichlet(LbBoundary):
"""Boundary condition for advection-diffusion problems that fixes the concentration at the obstacle.
Args:
concentration: value of the concentration which should be set.
name: optional name of the boundary.
"""
def __init__(self, concentration, name=None):
if name is None:
name = "Diffusion Dirichlet " + str(concentration)
self.concentration = concentration
super(DiffusionDirichlet, self).__init__(name)
def get_additional_code_nodes(self, lb_method):
"""Return a list of code nodes that will be added in the generated code before the index field loop.
Args:
lb_method: Lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`
Returns:
list containing LbmWeightInfo
"""
return [LbmWeightInfo(lb_method)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
w_dir = LbmWeightInfo.weight_of_direction(dir_symbol, lb_method)
return [Assignment(f_in(inv_dir[dir_symbol]),
2 * w_dir * self.concentration - f_out(dir_symbol))]
# end class DiffusionDirichlet
class NeumannByCopy(LbBoundary):
"""Neumann boundary condition which is implemented by coping the PDF values to achieve similar values at the fluid
and the boundary node"""
......@@ -482,14 +516,6 @@ class NeumannByCopy(LbBoundary):
return [Assignment(f_in(inv_dir[dir_symbol]), f_out(inv_dir[dir_symbol])),
Assignment(f_out[neighbour_offset](dir_symbol), f_out(dir_symbol))]
def __hash__(self):
# All boundaries of these class behave equal -> should also be equal
return hash("NeumannByCopy")
def __eq__(self, other):
return type(other) == NeumannByCopy
# end class NeumannByCopy
......@@ -504,7 +530,7 @@ class StreamInConstant(LbBoundary):
def __init__(self, constant, name=None):
super(StreamInConstant, self).__init__(name)
self._constant = constant
self.constant = constant
def get_additional_code_nodes(self, lb_method):
"""Return a list of code nodes that will be added in the generated code before the index field loop.
......@@ -519,13 +545,7 @@ class StreamInConstant(LbBoundary):
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
return [Assignment(f_in(inv_dir[dir_symbol]), self._constant),
Assignment(f_out[neighbour_offset](dir_symbol), self._constant)]
return [Assignment(f_in(inv_dir[dir_symbol]), self.constant),
Assignment(f_out[neighbour_offset](dir_symbol), self.constant)]
def __hash__(self):
# All boundaries of these class behave equal -> should also be equal
return hash("StreamInConstant")
def __eq__(self, other):
return type(other) == StreamInConstant
# end class StreamInConstant
lbmpy/continuous_distribution_measures.py
+ 39
− 7
View file @ d77a9064
......@@ -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
+ 9
− 6
View file @ d77a9064
......@@ -198,10 +198,11 @@ 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 FastCentralMomentTransform
from lbmpy.methods.momentbased.moment_transforms import PdfsToCentralMomentsByShiftMatrix
from lbmpy.methods.centeredcumulant.cumulant_transform import CentralMomentsToCumulantsByGeneratingFunc
from lbmpy.methods.creationfunctions import (
create_with_monomial_cumulants, create_with_polynomial_cumulants, create_with_default_polynomial_cumulants)
......@@ -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,))
......@@ -579,7 +582,7 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
'initial_velocity': None,
'galilean_correction': False, # only available for D3Q27 cumulant methods
'central_moment_transform_class': FastCentralMomentTransform,
'central_moment_transform_class': PdfsToCentralMomentsByShiftMatrix,
'cumulant_transform_class': CentralMomentsToCumulantsByGeneratingFunc,
'entropic': False,
......@@ -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
+ 1
− 1
View file @ d77a9064
......@@ -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
+ 35
− 43
View file @ d77a9064
......@@ -10,9 +10,13 @@ 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):
def get_weights(stencil, c_s_sq=sp.Rational(1, 3)):
q = len(stencil)
if c_s_sq != sp.Rational(1, 3) and c_s_sq != sp.Symbol("c_s") ** 2:
......@@ -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
+ 2
− 2
View file @ d77a9064
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
+ 2
− 2
View file @ d77a9064
......@@ -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 ----------------------------------------------------------
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment