diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index c207cada6018ec2b8bd497ff61a6a66e63637dab..558282733c85a6070944277f427f4e304e7ea063 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -233,6 +233,7 @@ class LBMConfig:
Adds the Cassons model according to https://doi.org/10.1007/s10955-005-8415-x
The parameters are set with the ``CassonsParameters`` dataclass.
"""
+
fluctuating: dict = False
"""
Enables fluctuating lattice Boltzmann by randomizing collision process.
@@ -244,6 +245,11 @@ class LBMConfig:
Temperature for fluctuating lattice Boltzmann methods.
"""
+ rr_parameters: Tuple = None
+ """
+ Parameters for the hybrid recursive collision operator sigma and velocity field
+ """
+
output: dict = field(default_factory=dict)
"""
A dictionary mapping macroscopic quantites e.g. the strings 'density' and 'velocity' to pystencils
@@ -728,7 +734,12 @@ def create_lb_method(lbm_config=None, **params):
assert len(relaxation_rates) >= 1, "Not enough relaxation rates"
method = create_srt(lbm_config.stencil, relaxation_rates[0], **common_params)
elif lbm_config.method == Method.RR:
- method = create_rr(lbm_config.stencil, relaxation_rates, **common_params)
+ sigma = sp.Integer(1)
+ velocity_field = None
+ if lbm_config.rr_parameters:
+ sigma = lbm_config.rr_parameters[0]
+ velocity_field = lbm_config.rr_parameters[1]
+ method = create_rr(lbm_config.stencil, relaxation_rates, sigma, velocity_field, **common_params)
elif lbm_config.method == Method.TRT:
assert len(relaxation_rates) >= 2, "Not enough relaxation rates"
method = create_trt(lbm_config.stencil, relaxation_rates[0], relaxation_rates[1], **common_params)
diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py
index 76b0c0a43b328010939c0ea150db756668bbf249..eecbffa9e6b12dda087b5b3abda74ced90c24072 100644
--- a/lbmpy/methods/creationfunctions.py
+++ b/lbmpy/methods/creationfunctions.py
@@ -250,13 +250,15 @@ def create_srt(stencil, relaxation_rate, continuous_equilibrium=True, **kwargs):
return create_with_discrete_maxwellian_equilibrium(stencil, rr_dict, **kwargs)
-def create_rr(stencil, relaxation_rates, **kwargs):
+def create_rr(stencil, relaxation_rates, sigma=sp.Integer(1), velocity_field=None, **kwargs):
r"""Creates a recursive regularised collision model
Args:
stencil: instance of :class:`lbmpy.stencils.LBStencil`
relaxation_rates: relaxation rate (inverse of the relaxation time)
usually called :math:`\omega` in LBM literature
+ sigma: switch between rr and hrr
+ velocity_field: velocity field to calculate the finite difference
Returns:
:class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
@@ -291,7 +293,8 @@ def create_rr(stencil, relaxation_rates, **kwargs):
order=equilibrium_order,
c_s_sq=c_s_sq)
- return RecursiveRegularisedLbMethod(stencil, equilibrium, rr_dict, cqc, force_model, zero_centered)
+ return RecursiveRegularisedLbMethod(stencil, equilibrium, rr_dict, cqc, force_model, zero_centered,
+ sigma, velocity_field)
def create_trt(stencil, relaxation_rate_even_moments, relaxation_rate_odd_moments,
diff --git a/lbmpy/methods/momentbased/recursiveregularisedmethod.py b/lbmpy/methods/momentbased/recursiveregularisedmethod.py
index 866de3614001d235353a8bb11ccb1ff0bfadc387..b48715b0ca8c4c0aea6f9b914d7fe5186c6cd1c1 100644
--- a/lbmpy/methods/momentbased/recursiveregularisedmethod.py
+++ b/lbmpy/methods/momentbased/recursiveregularisedmethod.py
@@ -53,7 +53,8 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
"""
def __init__(self, stencil, equilibrium, relaxation_dict,
- conserved_quantity_computation=None, force_model=None, zero_centered=False):
+ conserved_quantity_computation=None, force_model=None, zero_centered=False, sigma=sp.Integer(1),
+ velocity_field=None):
assert isinstance(conserved_quantity_computation, AbstractConservedQuantityComputation)
super(RecursiveRegularisedLbMethod, self).__init__(stencil)
@@ -63,6 +64,8 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
self._force_model = force_model
self._zero_centered = zero_centered
self._weights = None
+ self._sigma = sigma
+ self._velocity_field = velocity_field
@property
def force_model(self):
@@ -192,7 +195,8 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
additional_subexpressions=rr_sub_expressions,
include_force_terms=True,
conserved_quantity_equations=conserved_quantity_equations,
- pre_simplification=pre_simplification)
+ pre_simplification=pre_simplification,
+ sigma=self._sigma)
return ac
def set_zeroth_moment_relaxation_rate(self, relaxation_rate):
@@ -326,7 +330,8 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
additional_subexpressions: Iterable[Assignment] = None,
include_force_terms: bool = True,
conserved_quantity_equations: AssignmentCollection = None,
- pre_simplification: bool = False) -> LbmCollisionRule:
+ pre_simplification: bool = False,
+ sigma=sp.Integer(1)) -> LbmCollisionRule:
if additional_subexpressions is None:
additional_subexpressions = list()
@@ -351,6 +356,8 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
cs4 = c_s_sq * c_s_sq
u = sp.symbols("u_:3")
+ taup = rho * c_s_sq / omega1
+
RR = [sp.Symbol(f"M_{index.name[1:]}") for index in Ordering]
RReq = [sp.Symbol(f"eq_{index.name[1:]}") for index in Ordering]
# RMeq = [sp.Symbol(f"RMeq_{index.name[1:]}") for index in Ordering]
@@ -405,6 +412,27 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
M202 = Ordering.M202
M022 = Ordering.M022
+ gradient = []
+ sr = []
+
+ if self._velocity_field:
+ gradient = [Assignment(sp.Symbol("du_11"), self._velocity_field[1, 0, 0](0) - self._velocity_field[-1, 0, 0](0)),
+ Assignment(sp.Symbol("du_12"), self._velocity_field[0, 1, 0](0) - self._velocity_field[0, -1, 0](0)),
+ Assignment(sp.Symbol("du_13"), self._velocity_field[0, 0, 1](0) - self._velocity_field[0, 0, -1](0)),
+ Assignment(sp.Symbol("du_21"), self._velocity_field[1, 0, 0](1) - self._velocity_field[-1, 0, 0](1)),
+ Assignment(sp.Symbol("du_22"), self._velocity_field[0, 1, 0](1) - self._velocity_field[0, -1, 0](1)),
+ Assignment(sp.Symbol("du_23"), self._velocity_field[0, 0, 1](1) - self._velocity_field[0, 0, -1](1)),
+ Assignment(sp.Symbol("du_31"), self._velocity_field[1, 0, 0](2) - self._velocity_field[-1, 0, 0](2)),
+ Assignment(sp.Symbol("du_32"), self._velocity_field[0, 1, 0](2) - self._velocity_field[0, -1, 0](2)),
+ Assignment(sp.Symbol("du_33"), self._velocity_field[0, 0, 1](2) - self._velocity_field[0, 0, -1](2))]
+
+ sr = [Assignment(sp.Symbol("SR1"), 2.0 * sp.Symbol("du_11")),
+ Assignment(sp.Symbol("SR2"), 2.0 * sp.Symbol("du_22")),
+ Assignment(sp.Symbol("SR3"), 2.0 * sp.Symbol("du_33")),
+ Assignment(sp.Symbol("SR4"), sp.Symbol("du_12") + sp.Symbol("du_21")),
+ Assignment(sp.Symbol("SR5"), sp.Symbol("du_13") + sp.Symbol("du_31")),
+ Assignment(sp.Symbol("SR6"), sp.Symbol("du_23") + sp.Symbol("du_32"))]
+
moments = [Assignment(X_M1, f[FM00] + f[FMM0] + f[FMP0] + f[FM0M] + f[FM0P]),
Assignment(X_P1, f[FP00] + f[FPP0] + f[FPM0] + f[FP0P] + f[FP0M]),
Assignment(X_0,
@@ -449,12 +477,12 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
Assignment(RReq[M022], RReq[M020] * RReq[M002])]
neq = [
- Assignment(RRneq[M200], RR[M200] - RReq[M200]),
- Assignment(RRneq[M020], RR[M020] - RReq[M020]),
- Assignment(RRneq[M002], RR[M002] - RReq[M002]),
- Assignment(RRneq[M110], RR[M110] - RReq[M110]),
- Assignment(RRneq[M101], RR[M101] - RReq[M101]),
- Assignment(RRneq[M011], RR[M011] - RReq[M011]),
+ Assignment(RRneq[M200], (RR[M200] - RReq[M200]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR1"))),
+ Assignment(RRneq[M020], (RR[M020] - RReq[M020]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR2"))),
+ Assignment(RRneq[M002], (RR[M002] - RReq[M002]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR3"))),
+ Assignment(RRneq[M110], (RR[M110] - RReq[M110]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR4"))),
+ Assignment(RRneq[M101], (RR[M101] - RReq[M101]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR5"))),
+ Assignment(RRneq[M011], (RR[M011] - RReq[M011]) * sigma + (sp.Integer(1) - sigma) * (-taup * sp.Symbol("SR6"))),
Assignment(RRneq[M210], u[1] * RRneq[M200] + 2.0 * u[0] * RRneq[M110]),
Assignment(RRneq[M201], u[2] * RRneq[M200] + 2.0 * u[0] * RRneq[M101]),
@@ -572,7 +600,7 @@ class RecursiveRegularisedLbMethod(AbstractLbMethod):
Assignment(d[F0MM], sp.Rational(1, 2) * rho * (- RMcoll[M021] - RMcoll[M012]) + tmp6),
]
- subexpressions = list(additional_subexpressions) + cqe.all_assignments + moments + eq_moments + neq + col
+ subexpressions = list(additional_subexpressions) + gradient + sr + cqe.all_assignments + moments + eq_moments + neq + col
ac = AssignmentCollection(subexpressions=subexpressions, main_assignments=main)
ac = ac.new_without_unused_subexpressions()