Fluctuating LBM

- variances from temperature
- seed parameter
- block offset parameter

lbmpy/creationfunctions.py

@@ -75,9 +75,8 @@ LES methods:
 
 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
-
+- ``fluctuating``: enables fluctuating lattice Boltzmann by randomizing collision process.
+  Pass dictionary with parameters to  ``lbmpy.fluctuatinglb.add_fluctuations_to_collision_rule``
 
 
 Optimization Parameters
@@ -175,11 +174,12 @@ import sympy as sp
 
 import lbmpy.forcemodels as forcemodels
 from lbmpy.fieldaccess import (
-    AAEvenTimeStepAccessor, AAOddTimeStepAccessor, CollideOnlyInplaceAccessor, EsoTwistEvenTimeStepAccessor,
-    EsoTwistOddTimeStepAccessor, PdfFieldAccessor, PeriodicTwoFieldsAccessor, StreamPullTwoFieldsAccessor,
-    StreamPushTwoFieldsAccessor)
-from lbmpy.fluctuatinglb import fluctuation_correction, fluctuating_variance_equations
-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 add_fluctuations_to_collision_rule, 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
@@ -193,8 +193,6 @@ 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.rng import random_symbol
-from pystencils.simp.assignment_collection import SymbolGen
 from pystencils.stencil import have_same_entries
 
 
@@ -309,6 +307,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)]
@@ -323,25 +324,7 @@ def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
     collision_rule = simplification(collision_rule)
 
     if params['fluctuating']:
-        variances = params["fluctuating"]
-        if params["fluctuating"] is True:
-            variances = [v for v, _ in zip(iter(SymbolGen("variance")), lb_method.moments)]
-
-        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["temperature"] is not None:
-            temperature = sp.Symbol("temperature") if params["temperature"] is True else params["temperature"]
-            variance_equations = fluctuating_variance_equations(method=lb_method,
-                                                                temperature=temperature,
-                                                                c_s_sq=params["c_s_sq"],
-                                                                variances=variances)
-            collision_rule.subexpressions += variance_equations
-            collision_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
+        add_fluctuations_to_collision_rule(collision_rule, **params['fluctuating'])
 
     if params['entropic']:
         if params['smagorinsky']:

lbmpy/fieldaccess.py

@@ -9,6 +9,7 @@ from pystencils.stencil import inverse_direction
 
 __all__ = ['PdfFieldAccessor', 'CollideOnlyInplaceAccessor', 'StreamPullTwoFieldsAccessor',
            'AAEvenTimeStepAccessor', 'AAOddTimeStepAccessor',
+           'PeriodicTwoFieldsAccessor', 'StreamPushTwoFieldsAccessor',
            'EsoTwistEvenTimeStepAccessor', 'EsoTwistOddTimeStepAccessor',
            'visualize_pdf_field_accessor', 'visualize_field_mapping']
 

lbmpy/fluctuatinglb.py

@@ -2,41 +2,45 @@
 
 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 numpy as np
 import sympy as sp
 
 from lbmpy.moments import MOMENT_SYMBOLS
+from pystencils import Assignment, TypedSymbol
+from pystencils.rng import PhiloxFourFloats, random_symbol
 from pystencils.simp.assignment_collection import SymbolGen
-from pystencils import Assignment
 
 
-def fluctuating_variance_equations(method, temperature=sp.Symbol("fluctuating_temperature"),
-                                   c_s_sq=sp.Symbol("c_s") ** 2,
-                                   variances=SymbolGen("variance")):
+def add_fluctuations_to_collision_rule(collision_rule, temperature=None, variances=(),
+                                       block_offsets=(0, 0, 0), seed=TypedSymbol("seed", np.uint32),
+                                       rng_node=PhiloxFourFloats, c_s_sq=sp.Rational(1, 3)):
+    """"""
+    if not (temperature and not variances) or (temperature and variances):
+        raise ValueError("Fluctuating LBM: Pass either 'temperature' or 'variances'.")
+
+    method = collision_rule.method
+    if not variances:
+        variances = fluctuating_variance_from_temperature(method, temperature, c_s_sq)
+
+    rng_symbol_gen = random_symbol(collision_rule.subexpressions, seed,
+                                   rng_node=rng_node, dim=method.dim, offsets=block_offsets)
+    correction = fluctuation_correction(method, rng_symbol_gen, 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)
+
+
+def fluctuating_variance_from_temperature(method, temperature, c_s_sq=sp.Symbol("c_s") ** 2):
     """Produces variance equations according to (3.54) in Schiller08"""
     normalization_factors = abs(method.moment_matrix) * sp.Matrix(method.weights)
-    relaxation_rates_sqr = [r ** 2 for r in method.relaxation_rates]
     density = method.zeroth_order_equilibrium_moment_symbol
+    if method.conserved_quantity_computation.zero_centered_pdfs:
+        density += 1
     mu = temperature * density / c_s_sq
-
-    return [Assignment(v, sp.sqrt(mu * norm * (1 - rr)))
-            for v, norm, rr in zip(iter(variances), normalization_factors, relaxation_rates_sqr)]
+    return [sp.sqrt(mu * norm * (1 - (1 - rr) ** 2))
+            for norm, rr in zip(normalization_factors, method.relaxation_rates)]
 
 
 def method_with_rescaled_equilibrium_values(base_method):

lbmpy/methods/conservedquantitycomputation.py

@@ -243,10 +243,11 @@ def get_equations_for_zeroth_and_first_order_moment(stencil, symbolic_pdfs,
         \rho = \sum_{d \in S} f_d
         u_j = \sum_{d \in S} f_d u_jd
 
-    :param stencil: called :math:`S` above
-    :param symbolic_pdfs: called :math:`f` above
-    :param symbolic_zeroth_moment:  called :math:`\rho` above
-    :param symbolic_first_moments: called :math:`u` above
+    Args:
+        stencil: called :math:`S` above
+        symbolic_pdfs: called :math:`f` above
+        symbolic_zeroth_moment:  called :math:`\rho` above
+        symbolic_first_moments: called :math:`u` above
     """
     def filter_out_plus_terms(expr):
         result = 0

lbmpy_tests/test_fluctuation.py

+import numpy as np
+
 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=True,
-                        temperature=1e-9,
-                        kernel_params={'time_step': 1})
+                        fluctuating={'temperature': 1e-9},
+                        kernel_params={'time_step': 1, 'seed': 312},
+                        optimization={'cse_global': True})
 
     ch.run(10)
+    assert np.max(ch.velocity[:, :]) < 0.1