Support for fluctuating LB methods

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,

lbmpy/fluctuatinglb.py

+"""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

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):

lbmpy_tests/test_fluctuation.py

+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)