diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index acf9ee67b9fa924cc268c42dbed9adb95aea1f00..6c18b643d4e57b4106cb8bf0158262d4ef12b476 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -73,6 +73,12 @@ LES methods:
- ``smagorinsky=False``: set to Smagorinsky constant to activate turbulence model, ``omega_output_field`` can be set to
write out adapted relaxation rates
+Fluctuating LB:
+
+- ``fluctuating=(variance1, variance2, )``: enables fluctuating lattice Boltzmann by randomizing collision process.
+ Pass sequence of variances for each moment, or `True` to use symbolic variances
+
+
Optimization Parameters
^^^^^^^^^^^^^^^^^^^^^^^
@@ -169,11 +175,11 @@ import sympy as sp
import lbmpy.forcemodels as forcemodels
from lbmpy.fieldaccess import (
- AAEvenTimeStepAccessor, AAOddTimeStepAccessor, CollideOnlyInplaceAccessor,
- EsoTwistEvenTimeStepAccessor, EsoTwistOddTimeStepAccessor, PdfFieldAccessor,
- PeriodicTwoFieldsAccessor, StreamPullTwoFieldsAccessor, StreamPushTwoFieldsAccessor)
-from lbmpy.methods import (
- create_mrt3, create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc)
+ AAEvenTimeStepAccessor, AAOddTimeStepAccessor, CollideOnlyInplaceAccessor, EsoTwistEvenTimeStepAccessor,
+ EsoTwistOddTimeStepAccessor, PdfFieldAccessor, PeriodicTwoFieldsAccessor, StreamPullTwoFieldsAccessor,
+ StreamPushTwoFieldsAccessor)
+from lbmpy.fluctuatinglb import fluctuation_correction, method_with_rescaled_equilibrium_values
+from lbmpy.methods import create_mrt3, create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc
from lbmpy.methods.creationfunctions import create_generic_mrt
from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
from lbmpy.methods.entropic import add_entropy_condition, add_iterative_entropy_condition
@@ -187,7 +193,8 @@ from pystencils import Assignment, AssignmentCollection, create_kernel
from pystencils.cache import disk_cache_no_fallback
from pystencils.data_types import collate_types
from pystencils.field import Field, get_layout_of_array
-from pystencils.simp import add_subexpressions_for_field_reads
+from pystencils.rng import random_symbol
+from pystencils.simp.assignment_collection import SymbolGen
from pystencils.stencil import have_same_entries
@@ -302,6 +309,9 @@ def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
if rho_in is not None and isinstance(rho_in, Field):
rho_in = rho_in.center
+ if params['fluctuating']:
+ lb_method = method_with_rescaled_equilibrium_values(lb_method)
+
if u_in is not None:
density_rhs = sum(lb_method.pre_collision_pdf_symbols) if rho_in is None else rho_in
eqs = [Assignment(cqc.zeroth_order_moment_symbol, density_rhs)]
@@ -315,6 +325,14 @@ def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
collision_rule = simplification(collision_rule)
+ if params['fluctuating']:
+ variances = SymbolGen("variance") if params['fluctuating'] is True else params['fluctuating']
+ correction = fluctuation_correction(lb_method, random_symbol(collision_rule.subexpressions, dim=lb_method.dim),
+ variances)
+ for i, corr in enumerate(correction):
+ collision_rule.main_assignments[i] = Assignment(collision_rule.main_assignments[i].lhs,
+ collision_rule.main_assignments[i].rhs + corr)
+
if params['entropic']:
if params['smagorinsky']:
raise ValueError("Choose either entropic or smagorinsky")
@@ -509,6 +527,7 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
'entropic_newton_iterations': None,
'omega_output_field': None,
'smagorinsky': False,
+ 'fluctuating': False,
'output': {},
'velocity_input': None,
diff --git a/lbmpy/fluctuatinglb.py b/lbmpy/fluctuatinglb.py
new file mode 100644
index 0000000000000000000000000000000000000000..c0d4ab3b23fe7c75fa55226ec9f777f23ee9f634
--- /dev/null
+++ b/lbmpy/fluctuatinglb.py
@@ -0,0 +1,50 @@
+"""Functions for implementation of fluctuating (randomized) lattice Boltzmann
+
+to generate a fluctuating LBM the equilibrium moment values have to be scaled and an additive (random)
+correction term is added to the collision rule
+
+
+Usage example:
+
+ >>> from lbmpy.session import *
+ >>> from pystencils.rng import random_symbol
+ >>> method = create_lb_method(stencil='D2Q9', method='MRT')
+ >>> rescaled_method = method_with_rescaled_equilibrium_values(method)
+ >>> cr = create_lb_collision_rule(lb_method=rescaled_method)
+ >>> correction = fluctuation_correction(rescaled_method,
+ ... rng_generator=random_symbol(cr.subexpressions, dim=method.dim))
+ >>> for i, corr in enumerate(correction):
+ ... cr.main_assignments[i] = ps.Assignment(cr.main_assignments[i].lhs,
+ ... cr.main_assignments[i].rhs + corr)
+ >>>
+
+"""
+import sympy as sp
+
+from lbmpy.moments import MOMENT_SYMBOLS
+from pystencils.simp.assignment_collection import SymbolGen
+
+
+def method_with_rescaled_equilibrium_values(base_method):
+ """Re-scales the equilibrium moments by 1 / sqrt(M*w) with moment matrix M and weights w"""
+ from lbmpy.creationfunctions import create_lb_method_from_existing
+
+ sig_k = abs(base_method.moment_matrix) * sp.Matrix(base_method.weights)
+
+ def modification_rule(moment, eq, rr):
+ i = base_method.moments.index(moment)
+ return moment, eq / sp.sqrt(sig_k[i]), rr
+
+ return create_lb_method_from_existing(base_method, modification_rule)
+
+
+def fluctuation_correction(method, rng_generator, variances=SymbolGen("variance")):
+ """Returns additive correction terms to be added to the the collided pdfs"""
+ conserved_moments = {sp.sympify(1), *MOMENT_SYMBOLS}
+
+ # A diagonal matrix containing the random fluctuations
+ random_matrix = sp.Matrix([0 if m in conserved_moments else next(rng_generator) for m in method.moments])
+ random_variance = sp.diag(*[v for v, _ in zip(iter(variances), method.moments)])
+
+ # corrections are applied in real space hence we need to convert to real space here
+ return method.moment_matrix.inv() * random_variance * random_matrix
diff --git a/lbmpy/forcemodels.py b/lbmpy/forcemodels.py
index 7b116e462507e5372e02675b414747f2642cedbb..017927c37c8cb27f962497ab24770f35a11329de 100644
--- a/lbmpy/forcemodels.py
+++ b/lbmpy/forcemodels.py
@@ -91,7 +91,7 @@ import sympy as sp
from lbmpy.relaxationrates import get_shear_relaxation_rate
-class Simple(object):
+class Simple:
r"""
A simple force model which introduces the following additional force term in the
collision process :math:`\frac{w_q}{c_s^2} c_{qi} \; a_i` (often: force = rho * acceleration)
@@ -114,7 +114,7 @@ class Simple(object):
return default_velocity_shift(density, self._force)
-class Luo(object):
+class Luo:
r"""Force model by Luo :cite:`luo1993lattice`.
Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity.
@@ -136,7 +136,7 @@ class Luo(object):
return default_velocity_shift(density, self._force)
-class Guo(object):
+class Guo:
r"""
Force model by Guo :cite:`guo2002discrete`
Adapts the calculation of the macroscopic velocity as well as the equilibrium velocity (both shifted by F/2)!
@@ -158,7 +158,7 @@ class Guo(object):
return default_velocity_shift(density, self._force)
-class Buick(object):
+class Buick:
r"""
This force model :cite:`buick2000gravity` has a force term with zero second moment.
It is suited for incompressible lattice models.
@@ -181,7 +181,7 @@ class Buick(object):
return default_velocity_shift(density, self._force)
-class EDM(object):
+class EDM:
r"""Exact differencing force model"""
def __init__(self, force):
diff --git a/lbmpy_tests/test_fluctuation.py b/lbmpy_tests/test_fluctuation.py
new file mode 100644
index 0000000000000000000000000000000000000000..f3e45f5b20cb8a3aeda5646edbde5dda1bd9764c
--- /dev/null
+++ b/lbmpy_tests/test_fluctuation.py
@@ -0,0 +1,8 @@
+from lbmpy.scenarios import create_channel
+
+def test_fluctuating_generation_pipeline():
+ ch = create_channel((40, 10), method='mrt3', relaxation_rates=[1.5, 1, 1], force=1e-5,
+ fluctuating=[1e-3] * 9,
+ kernel_params={'time_step': 1})
+
+ ch.run(10)