fluctuating MRT

lbmpy/creationfunctions.py

@@ -175,7 +175,7 @@ from pystencils.field import get_layout_of_array, Field
 from pystencils import create_kernel, Assignment
 from lbmpy.turbulence_models import add_smagorinsky_model
 from lbmpy.methods import create_srt, create_trt, create_mrt_orthogonal, create_trt_kbc, \
-    create_mrt_raw, create_mrt3
+    create_mrt_raw, create_mrt3, create_fluctuating_mrt
 from lbmpy.methods.entropic import add_iterative_entropy_condition, add_entropy_condition
 from lbmpy.methods.entropic_eq_srt import create_srt_entropic
 from lbmpy.relaxationrates import relaxation_rate_from_magic_number
@@ -388,8 +388,24 @@ def create_lb_method(**params):
         method = create_mrt_orthogonal(stencil_entries, relaxation_rate_getter, **common_params)
     elif method_name.lower() == 'mrt_raw':
         method = create_mrt_raw(stencil_entries, relaxation_rates, **common_params)
+    elif method_name.lower() == 'fluctuating_mrt':
+        next_relaxation_rate = [0]
+
+        def relaxation_rate_getter(_):
+            res = relaxation_rates[next_relaxation_rate[0]]
+            next_relaxation_rate[0] += 1
+            return res
+        method1 = create_mrt_orthogonal(stencil_entries, relaxation_rate_getter, **common_params)
+        if not 'random_number_symbols' in params:
+            params['random_number_symbols'] = sp.symbols("rand_:15")
+        method = create_fluctuating_mrt(method1, common_params['maxwellian_moments'], params['random_number_symbols'])
     elif method_name.lower() == 'mrt3':
         method = create_mrt3(stencil_entries, relaxation_rates, **common_params)
+    elif method_name.lower() == 'fluctuating_mrt3':
+        method1 = create_mrt3(stencil_entries, relaxation_rates, **common_params)
+        if not 'random_number_symbols' in params:
+            params['random_number_symbols'] = sp.symbols("rand_:15")
+        method = create_fluctuating_mrt(method1, common_params['maxwellian_moments'], params['random_number_symbols'])
     elif method_name.lower().startswith('trt-kbc-n'):
         if have_same_entries(stencil_entries, get_stencil("D2Q9")):
             dim = 2

lbmpy/methods/__init__.py

@@ -3,10 +3,11 @@ from lbmpy.methods.momentbased import RelaxationInfo, AbstractLbMethod, Abstract
 from .conservedquantitycomputation import DensityVelocityComputation
 from lbmpy.methods.creationfunctions import create_srt, create_trt, create_trt_with_magic_number, create_trt_kbc, \
     create_mrt_orthogonal, create_mrt_raw, create_mrt3, \
+    create_fluctuating_mrt,\
     create_with_continuous_maxwellian_eq_moments, create_with_discrete_maxwellian_eq_moments
 
 __all__ = ['RelaxationInfo', 'AbstractLbMethod',
            'AbstractConservedQuantityComputation', 'DensityVelocityComputation', 'MomentBasedLbMethod',
            'create_srt', 'create_trt', 'create_trt_with_magic_number', 'create_trt_kbc',
-           'create_mrt_orthogonal', 'create_mrt_raw', 'create_mrt3',
+           'create_mrt_orthogonal', 'create_mrt_raw', 'create_mrt3', 'create_fluctuating_mrt',
            'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments']

lbmpy/methods/creationfunctions.py

@@ -219,6 +219,81 @@ def create_trt_with_magic_number(stencil, relaxation_rate, magic_number=sp.Ratio
                       relaxation_rate_odd_moments=rr_odd, **kwargs)
 
 
+def create_fluctuating_mrt(base_mrt_method, maxwellian_moments=False, random_number_symbols=sp.symbols("rand_:15")):
+    """
+    Creates a fluctuating mrt from a normal mrt method passed in by adding a fluctuating force term.
+    :param base_mrt_method:
+    :return:
+    """
+    class FluctuatingForce(object):
+        def __init__(self, random_numbers, random_variances):
+            self.random_numbers = random_numbers
+            self.random_variances = random_variances
+
+        def __call__(self, lb_method, **kwargs):
+            # Random fluctuations on all but the 4 conserved modes
+            assert (len(self.random_numbers) == (len(lb_method.stencil) - 4))
+
+            moment_matrix = lb_method.moment_matrix
+
+            # A diagonal matrix containing the random fluctuations
+            random_matrix = sp.Matrix([0, 0, 0, 0, *self.random_numbers])
+            random_variance = sp.diag(*self.random_variances)
+
+            # Forces are applied in real space hence we need to convert to real space here
+            return moment_matrix.inv() * random_variance * random_matrix
+
+    mu = sp.Symbol("mu", positive=True)
+    random_variances = sp.symbols("phi_:19")
+    variance_eqs = []
+
+    for i in range(0, 4):
+        variance_eqs.append(sp.Eq(random_variances[i], 0))
+
+    for i in range(4, 19):
+        r = base_mrt_method.relaxation_rates[i]
+        variance_eqs.append(sp.Eq(random_variances[i], sp.sqrt(mu) * sp.sqrt(2 * r - r * r)))
+
+    # Weighted length of rows
+    sigk = (abs(base_mrt_method.moment_matrix) * sp.Matrix(base_mrt_method.weights))
+
+    inv_sqrt_sigk = sigk.applyfunc(lambda x: (1 / sp.sqrt(x)))
+
+    mm_to_rr_dict = OrderedDict()
+    i = 0
+    for moment, value in base_mrt_method.relaxation_info_dict.items():
+        modified_moment = moment * inv_sqrt_sigk[i]
+        mm_to_rr_dict[modified_moment] = value.relaxation_rate
+        i += 1
+
+    random_variances = sp.symbols("phi_:19")
+
+    force_model = FluctuatingForce(random_number_symbols, random_variances)
+    if maxwellian_moments:
+        fluctuating_base = create_with_continuous_maxwellian_eq_moments(base_mrt_method.stencil, mm_to_rr_dict,
+                                                                        force_model=force_model)
+    else:
+        fluctuating_base = create_with_discrete_maxwellian_eq_moments(base_mrt_method.stencil, mm_to_rr_dict,
+                                                                      force_model=force_model)
+
+    """
+    TODO: Add new relaxation rate equations
+    
+    nu = sp.Symbol("nu")
+    nub = sp.Symbol("nu_b")
+    relaxation_rates = fluctuating_base.relaxation_rates
+    relaxation_equations = [sp.Eq(sp.Symbol("omega_0"), 0)]
+    relaxation_equations += [sp.Eq(relaxation_rates[4], -1*((2*nu-1)/(2*nu+1) - 1))]
+    relaxation_equations += [sp.Eq(relaxation_rates[9], -1*((3*nub-1)/(3*nub+1) -1))]
+    relaxation_equations += [sp.Eq(sp.Symbol("omega_2"), 1)]
+    relaxation_equations += [sp.Eq(sp.Symbol("omega_3"), 1)]
+    relaxation_equations += [sp.Eq(sp.Symbol("omega_5"), 1)]
+    relaxation_equations += [sp.Eq(sp.Symbol("omega_6"), 1)]
+    """
+
+    return fluctuating_base
+
+
 def create_mrt_raw(stencil, relaxation_rates, maxwellian_moments=False, **kwargs):
     """Creates a MRT method using non-orthogonalized moments."""
     if isinstance(stencil, str):