diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 4ed8d633dd7f749af7fb8ecf7caf18bdbeb6c80d..f708e7f98876b048d8508b7afaff8934929356a1 100644
--- a/lbmpy/creationfunctions.py
+++ b/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
diff --git a/lbmpy/methods/__init__.py b/lbmpy/methods/__init__.py
index 3e068b6f0eef5c41de27d172f4cb3cc3c5ee1135..b6552a476c39c1528a3f757a8331ce96f44e1e0e 100644
--- a/lbmpy/methods/__init__.py
+++ b/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']
diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py
index e910fc413b5b37e4dbe99c3368055b60987ae6d5..e835fe1db6c7a82a3e00a2457d77476c3b4cac09 100644
--- a/lbmpy/methods/creationfunctions.py
+++ b/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):