diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py index 83a1ecfcc44777e17a59f88133ce2ce590516908..73769bf97feaf33479cdb69732c371765c89babe 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 b43300d14faaf1747d3a77147b69383488251b83..3eee11f027f1d01db1f796a7390daf2162a89bcf 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 0f018a649a853cc2c42cacb41474b882bde8799e..650ece1359046ad054aa50e8559f8cb0ff530e9f 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 be0e953f22b7b221bc4fbeb0bf4bcafd2fab7351..7ac94f42b5cf649147cbefa5830d2221c628cc14 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])