From adefba32e0b906416db7514df529f8f09efbc408 Mon Sep 17 00:00:00 2001
From: Michael Kuron <mkuron@icp.uni-stuttgart.de>
Date: Mon, 18 Nov 2019 16:38:32 +0100
Subject: [PATCH] Fluctuating MRT: requires weighted-orthogonal method
---
lbmpy/fluctuatinglb.py | 3 +++
lbmpy/methods/momentbased.py | 11 +++++++++++
lbmpy_tests/test_fluctuation.py | 7 +++----
lbmpy_tests/test_momentbased_methods_equilibrium.py | 7 +++++--
4 files changed, 22 insertions(+), 6 deletions(-)
diff --git a/lbmpy/fluctuatinglb.py b/lbmpy/fluctuatinglb.py
index 86aebd48..069fa999 100644
--- a/lbmpy/fluctuatinglb.py
+++ b/lbmpy/fluctuatinglb.py
@@ -25,6 +25,9 @@ def add_fluctuations_to_collision_rule(collision_rule, temperature=None, amplitu
if block_offsets == 'walberla':
block_offsets = tuple(TypedSymbol("block_offset_{}".format(i), np.uint32) for i in range(3))
+ if not method.is_weighted_orthogonal:
+ raise ValueError("Fluctuations can only be added to weighted-orthogonal methods")
+
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, amplitudes)
diff --git a/lbmpy/methods/momentbased.py b/lbmpy/methods/momentbased.py
index 67d0ec2d..50f0ef87 100644
--- a/lbmpy/methods/momentbased.py
+++ b/lbmpy/methods/momentbased.py
@@ -2,6 +2,7 @@ from collections import OrderedDict
import sympy as sp
+from lbmpy.maxwellian_equilibrium import get_weights
from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
from lbmpy.moments import MOMENT_SYMBOLS, moment_matrix
@@ -124,6 +125,16 @@ class MomentBasedLbMethod(AbstractLbMethod):
def moment_matrix(self):
return moment_matrix(self.moments, self.stencil)
+ @property
+ def is_orthogonal(self):
+ return (self.moment_matrix * self.moment_matrix.T).is_diagonal()
+
+ @property
+ def is_weighted_orthogonal(self):
+ w = get_weights(self.stencil, sp.Rational(1, 3))
+ return (sp.matrix_multiply_elementwise(self.moment_matrix, sp.Matrix([w] * len(w))) * self.moment_matrix.T
+ ).is_diagonal()
+
def __getstate__(self):
# Workaround for a bug in joblib
self._momentToRelaxationInfoDictToPickle = [i for i in self._momentToRelaxationInfoDict.items()]
diff --git a/lbmpy_tests/test_fluctuation.py b/lbmpy_tests/test_fluctuation.py
index be886c28..18e4f358 100644
--- a/lbmpy_tests/test_fluctuation.py
+++ b/lbmpy_tests/test_fluctuation.py
@@ -4,10 +4,9 @@ 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={'temperature': 1e-9},
- kernel_params={'time_step': 1, 'seed': 312},
+ ch = create_channel((10, 10, 10), stencil='D3Q19', method='mrt', relaxation_rates=[1.5] * 7, force=1e-5,
+ 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
+ assert np.max(ch.velocity[:, :, :]) < 0.1
diff --git a/lbmpy_tests/test_momentbased_methods_equilibrium.py b/lbmpy_tests/test_momentbased_methods_equilibrium.py
index 843685bc..f2698f3c 100644
--- a/lbmpy_tests/test_momentbased_methods_equilibrium.py
+++ b/lbmpy_tests/test_momentbased_methods_equilibrium.py
@@ -70,7 +70,10 @@ def test_relaxation_rate_setter():
def test_mrt_orthogonal():
m = create_mrt_orthogonal(get_stencil("D2Q9"), maxwellian_moments=True)
- assert (m.moment_matrix * m.moment_matrix.T).is_diagonal()
+ assert m.is_orthogonal
+
+ m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True)
+ assert m.is_weighted_orthogonal
m = create_mrt_orthogonal(get_stencil("D3Q27"), maxwellian_moments=True)
- assert (m.moment_matrix * m.moment_matrix.T).is_diagonal()
+ assert m.is_orthogonal
--
GitLab