From 1175fb4a44144aebee80339b71c21837b9f51ca1 Mon Sep 17 00:00:00 2001
From: Markus Holzer <markus.holzer@fau.de>
Date: Thu, 15 Jun 2023 15:37:23 +0200
Subject: [PATCH] Introduce Central moment forcing again
---
lbmpy/creationfunctions.py | 8 ++++++--
lbmpy/enums.py | 4 ++++
lbmpy/forcemodels.py | 22 ++++++++++++++++++++++
lbmpy_tests/test_force.py | 3 +--
4 files changed, 33 insertions(+), 4 deletions(-)
diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 83a1ecfc..73769bf9 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -412,7 +412,10 @@ class LBMConfig:
force_not_zero = True
if self.force_model is None and force_not_zero:
- self.force_model = forcemodels.Guo(self.force[:self.stencil.D])
+ if self.method == Method.CUMULANT:
+ self.force_model = forcemodels.CentralMoment(self.force[:self.stencil.D])
+ else:
+ self.force_model = forcemodels.Guo(self.force[:self.stencil.D])
force_model_dict = {
'simple': forcemodels.Simple,
@@ -424,7 +427,8 @@ class LBMConfig:
'edm': forcemodels.EDM,
'kupershtokh': forcemodels.EDM,
'he': forcemodels.He,
- 'shanchen': forcemodels.ShanChen
+ 'shanchen': forcemodels.ShanChen,
+ 'centralmoment': forcemodels.CentralMoment
}
if isinstance(self.force_model, str):
diff --git a/lbmpy/enums.py b/lbmpy/enums.py
index b43300d1..3eee11f0 100644
--- a/lbmpy/enums.py
+++ b/lbmpy/enums.py
@@ -217,3 +217,7 @@ class ForceModel(Enum):
"""
See :class:`lbmpy.forcemodels.ShanChen`
"""
+ CENTRALMOMENT = auto()
+ """
+ See :class:`lbmpy.forcemodels`
+ """
diff --git a/lbmpy/forcemodels.py b/lbmpy/forcemodels.py
index 0f018a64..650ece13 100644
--- a/lbmpy/forcemodels.py
+++ b/lbmpy/forcemodels.py
@@ -210,6 +210,28 @@ class Simple(AbstractForceModel):
return before, after
+class CentralMoment(AbstractForceModel):
+ r"""
+ A force model that only shifts the macroscopic and equilibrium velocities and does not introduce a forcing term to
+ the collision process. Forcing is then applied through relaxation of the first central moments in the shifted
+ frame of reference (cf. https://doi.org/10.1016/j.camwa.2015.05.001).
+ """
+ def __call__(self, lb_method):
+ raise ValueError("This force model can only be combined with the Cumulant collision model")
+
+ def symmetric_central_moment_forcing(self, lb_method, *_):
+ before = sp.Matrix([0, ] * lb_method.stencil.Q)
+ after = sp.Matrix([0, ] * lb_method.stencil.Q)
+ for i, sf in enumerate(self.symbolic_force_vector):
+ before[i + 1] = sp.Rational(1, 2) * sf
+ after[i + 1] = sp.Rational(1, 2) * sf
+
+ return before, after
+
+ def equilibrium_velocity_shift(self, density):
+ return default_velocity_shift(density, self.symbolic_force_vector)
+
+
class Luo(AbstractForceModel):
r"""Force model by Luo :cite:`luo1993lattice`.
diff --git a/lbmpy_tests/test_force.py b/lbmpy_tests/test_force.py
index be0e953f..7ac94f42 100644
--- a/lbmpy_tests/test_force.py
+++ b/lbmpy_tests/test_force.py
@@ -12,10 +12,9 @@ from lbmpy.macroscopic_value_kernels import macroscopic_values_setter, macroscop
from lbmpy.moments import (is_bulk_moment, moments_up_to_component_order,
exponents_to_polynomial_representations, exponent_tuple_sort_key)
from lbmpy.stencils import LBStencil
-from lbmpy.updatekernels import create_stream_pull_with_output_kernel
# all force models available are defined in the ForceModel enum, but Cumulant is not a "real" force model
-force_models = [f for f in ForceModel]
+force_models = [f for f in ForceModel if f is not ForceModel.CENTRALMOMENT]
@pytest.mark.parametrize("method_enum", [Method.SRT, Method.TRT, Method.MRT, Method.CUMULANT])
--
GitLab