diff --git a/conftest.py b/conftest.py
index e71dff56416d4fa2a5c01eab39e1853b90fe1969..ebd72bdc77cc352960a98b4a763feb18ce76628a 100644
--- a/conftest.py
+++ b/conftest.py
@@ -60,7 +60,10 @@ try:
except ImportError:
collect_ignore += [os.path.join(SCRIPT_FOLDER, "lbmpy_tests/benchmark"),
os.path.join(SCRIPT_FOLDER,
- "lbmpy_tests/full_scenarios/kida_vortex_flow/scenario_kida_vortex_flow.py")]
+ "lbmpy_tests/full_scenarios/kida_vortex_flow/scenario_kida_vortex_flow.py"),
+ os.path.join(SCRIPT_FOLDER, "lbmpy_tests/full_scenarios/shear_wave/scenario_shear_wave.py"),
+ os.path.join(SCRIPT_FOLDER, "lbmpy_tests/test_json_serializer.py"),
+ os.path.join(SCRIPT_FOLDER, "lbmpy/db.py")]
if platform.system().lower() == 'windows':
collect_ignore += [os.path.join(SCRIPT_FOLDER, "lbmpy_tests/test_quicktests.py")]
diff --git a/doc/sphinx/lbmpy.bib b/doc/sphinx/lbmpy.bib
index 88198b2dcd4685c7024604b0b878ab0cd95b94c0..896a6a2a8deef32d432b8d004b43938c1c19fa10 100644
--- a/doc/sphinx/lbmpy.bib
+++ b/doc/sphinx/lbmpy.bib
@@ -110,6 +110,24 @@ issn = {0898-1221},
doi = {10.1016/j.camwa.2015.05.001},
}
+@article{geier2017,
+author = {Geier, Martin and Pasquali, Andrea and Sch{\"{o}}nherr, Martin},
+title = {Parametrization of the Cumulant Lattice Boltzmann Method for Fourth Order Accurate Diffusion Part I},
+year = {2017},
+issue_date = {November 2017},
+publisher = {Academic Press Professional, Inc.},
+address = {USA},
+volume = {348},
+number = {C},
+issn = {0021-9991},
+url = {https://doi.org/10.1016/j.jcp.2017.05.040},
+doi = {10.1016/j.jcp.2017.05.040},
+journal = {J. Comput. Phys.},
+month = {nov},
+pages = {862–888},
+numpages = {27}
+}
+
@Article{Coreixas2019,
title = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
author = {Coreixas, Christophe and Chopard, Bastien and Latt, Jonas},
diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 73769bf97feaf33479cdb69732c371765c89babe..e27501e2b501a1f3097fc2154dd962da66bedc4e 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -87,7 +87,7 @@ from pystencils.simp import sympy_cse, SimplificationStrategy
from lbmpy.methods.abstractlbmethod import LbmCollisionRule, AbstractLbMethod
from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
-# Filter out JobLib warnings. They are not usefull for use:
+# Filter out JobLib warnings. They are not useful for use:
# https://github.com/joblib/joblib/issues/683
filterwarnings("ignore", message="Persisting input arguments took")
@@ -99,7 +99,7 @@ class LBMConfig:
"""
stencil: lbmpy.stencils.LBStencil = LBStencil(Stencil.D2Q9)
"""
- All stencils are defined in :class:`lbmpy.enums.Stencil`. From that :class:`lbmpy.stencils.LBStenil`
+ All stencils are defined in :class:`lbmpy.enums.Stencil`. From that :class:`lbmpy.stencils.LBStencil`
class will be created
"""
method: Method = Method.SRT
@@ -197,6 +197,13 @@ class LBMConfig:
Special correction for D3Q27 cumulant LBMs. For Details see
:mod:`lbmpy.methods.cumulantbased.galilean_correction`
"""
+ fourth_order_correction: Union[float, bool] = False
+ """
+ Special correction for rendering D3Q27 cumulant LBMs fourth-order accurate in diffusion. For Details see
+ :mod:`lbmpy.methods.cumulantbased.fourth_order_correction`. If set to `True`, the fourth-order correction is
+ employed without limiters (or more precisely with a very high limiter, practically disabling the limiters). If this
+ variable is set to a number, the latter is used for the limiters (uniformly for omega_3, omega_4 and omega_5).
+ """
collision_space_info: CollisionSpaceInfo = None
"""
Information about the LB method's collision space (see :class:`lbmpy.methods.creationfunctions.CollisionSpaceInfo`)
@@ -646,6 +653,19 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
from lbmpy.methods.cumulantbased import add_galilean_correction
collision_rule = add_galilean_correction(collision_rule)
+ if lbm_config.fourth_order_correction:
+ from lbmpy.methods.cumulantbased import add_fourth_order_correction
+
+ # must provide a second relaxation rate in implementation; defaults to 1
+ if len(lbm_config.relaxation_rates) == 1:
+ lbm_config.relaxation_rates.append(1)
+
+ cumulant_limiter = 1e6 if lbm_config.fourth_order_correction is True else lbm_config.fourth_order_correction
+ collision_rule = add_fourth_order_correction(collision_rule=collision_rule,
+ shear_relaxation_rate=lbm_config.relaxation_rates[0],
+ bulk_relaxation_rate=lbm_config.relaxation_rates[1],
+ limiter=cumulant_limiter)
+
if lbm_config.entropic:
if lbm_config.smagorinsky or lbm_config.cassons:
raise ValueError("Choose either entropic, smagorinsky or cassons")
@@ -752,11 +772,25 @@ def create_lb_method(lbm_config=None, **params):
method_nr = lbm_config.method.name[-1]
method = create_trt_kbc(dim, relaxation_rates[0], relaxation_rates[1], 'KBC-N' + method_nr, **common_params)
elif lbm_config.method == Method.CUMULANT:
+ if lbm_config.fourth_order_correction:
+ if lbm_config.stencil.D != 3 and lbm_config.stencil.Q != 27:
+ raise ValueError("Fourth-order correction can only be applied to D3Q27 cumulant methods.")
+
+ assert len(relaxation_rates) <= 2, "Optimal parametrisation for fourth-order cumulants needs either one " \
+ "or two relaxation rates, associated with the shear (and bulk) " \
+ "viscosity. All other relaxation rates are automatically chosen " \
+ "optimally"
+
+ # define method in terms of symbolic relaxation rates and assign optimal values later
+ from lbmpy.methods.cumulantbased.fourth_order_correction import FOURTH_ORDER_RELAXATION_RATE_SYMBOLS
+ relaxation_rates = FOURTH_ORDER_RELAXATION_RATE_SYMBOLS
+
if lbm_config.nested_moments is not None:
method = create_cumulant(
lbm_config.stencil, relaxation_rates, lbm_config.nested_moments, **cumulant_params)
else:
method = create_with_default_polynomial_cumulants(lbm_config.stencil, relaxation_rates, **cumulant_params)
+
elif lbm_config.method == Method.MONOMIAL_CUMULANT:
method = create_with_monomial_cumulants(lbm_config.stencil, relaxation_rates, **cumulant_params)
else:
diff --git a/lbmpy/db.py b/lbmpy/db.py
new file mode 100644
index 0000000000000000000000000000000000000000..bca65fb182a125600b2e309dd2d5784f6d4d8c47
--- /dev/null
+++ b/lbmpy/db.py
@@ -0,0 +1,47 @@
+import json
+import six
+import inspect
+
+from pystencils.runhelper.db import PystencilsJsonEncoder
+from pystencils.simp import SimplificationStrategy
+from lbmpy import LBStencil, Method, CollisionSpace
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation
+from lbmpy.methods import CollisionSpaceInfo
+from lbmpy.forcemodels import AbstractForceModel
+from lbmpy.non_newtonian_models import CassonsParameters
+
+
+class LbmpyJsonEncoder(PystencilsJsonEncoder):
+
+ def default(self, obj):
+ if isinstance(obj, (LBMConfig, LBMOptimisation, CollisionSpaceInfo, CassonsParameters)):
+ return obj.__dict__
+ if isinstance(obj, (LBStencil, Method, CollisionSpace)):
+ return obj.name
+ if isinstance(obj, AbstractForceModel):
+ return obj.__class__.__name__
+ if isinstance(obj, SimplificationStrategy):
+ return obj.__str__()
+ if inspect.isclass(obj):
+ return obj.__name__
+ return PystencilsJsonEncoder.default(self, obj)
+
+
+class LbmpyJsonSerializer(object):
+
+ @classmethod
+ def serialize(cls, data):
+ if six.PY3:
+ if isinstance(data, bytes):
+ return json.dumps(data.decode('utf-8'), cls=LbmpyJsonEncoder, ensure_ascii=False).encode('utf-8')
+ else:
+ return json.dumps(data, cls=LbmpyJsonEncoder, ensure_ascii=False).encode('utf-8')
+ else:
+ return json.dumps(data, cls=LbmpyJsonEncoder, ensure_ascii=False).encode('utf-8')
+
+ @classmethod
+ def deserialize(cls, data):
+ if six.PY3:
+ return json.loads(data.decode('utf-8'))
+ else:
+ return json.loads(data.decode('utf-8'))
diff --git a/lbmpy/methods/abstractlbmethod.py b/lbmpy/methods/abstractlbmethod.py
index 03747ac9a0808f545982490e35cb34dee1852394..096a3dd6f264a0a15b12b813a24a251d2abd1514 100644
--- a/lbmpy/methods/abstractlbmethod.py
+++ b/lbmpy/methods/abstractlbmethod.py
@@ -67,8 +67,8 @@ class AbstractLbMethod(abc.ABC):
return d
@property
- def subs_dict_relxation_rate(self):
- """returns a dictonary which maps the replaced numerical relaxation rates to its original numerical value"""
+ def subs_dict_relaxation_rate(self):
+ """returns a dictionary which maps the replaced numerical relaxation rates to its original numerical value"""
result = dict()
for i in range(self._stencil.Q):
result[self.symbolic_relaxation_matrix[i, i]] = self.relaxation_matrix[i, i]
diff --git a/lbmpy/methods/cumulantbased/__init__.py b/lbmpy/methods/cumulantbased/__init__.py
index 939617543159a0d10584e170acfb9030f72a16bb..604b48038a5e75a44410fab1c3c6090c574d900f 100644
--- a/lbmpy/methods/cumulantbased/__init__.py
+++ b/lbmpy/methods/cumulantbased/__init__.py
@@ -1,4 +1,5 @@
from .cumulantbasedmethod import CumulantBasedLbMethod
from .galilean_correction import add_galilean_correction
+from .fourth_order_correction import add_fourth_order_correction
-__all__ = ['CumulantBasedLbMethod', 'add_galilean_correction']
+__all__ = ['CumulantBasedLbMethod', 'add_galilean_correction', 'add_fourth_order_correction']
diff --git a/lbmpy/methods/cumulantbased/cumulantbasedmethod.py b/lbmpy/methods/cumulantbased/cumulantbasedmethod.py
index c0bd07008daf2afb4e0a6867c54544512602dd0a..6d013958890056be5704544700cdabe6c3250de6 100644
--- a/lbmpy/methods/cumulantbased/cumulantbasedmethod.py
+++ b/lbmpy/methods/cumulantbased/cumulantbasedmethod.py
@@ -23,7 +23,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
This class implements cumulant-based lattice boltzmann methods which relax all the non-conserved quantities
as either monomial or polynomial cumulants. It is mostly inspired by the work presented in :cite:`geier2015`.
- This method is implemented modularily as the transformation from populations to central moments to cumulants
+ This method is implemented modularly as the transformation from populations to central moments to cumulants
is governed by subclasses of :class:`lbmpy.moment_transforms.AbstractMomentTransform`
which can be specified by constructor argument. This allows the selection of the most efficient transformation
for a given setup.
@@ -124,7 +124,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
@property
def zeroth_order_equilibrium_moment_symbol(self, ):
"""Returns a symbol referring to the zeroth-order moment of this method's equilibrium distribution,
- which is the area under it's curve
+ which is the area under its curve
(see :attr:`lbmpy.equilibrium.AbstractEquilibrium.zeroth_order_moment_symbol`)."""
return self._equilibrium.zeroth_order_moment_symbol
diff --git a/lbmpy/methods/cumulantbased/fourth_order_correction.py b/lbmpy/methods/cumulantbased/fourth_order_correction.py
new file mode 100644
index 0000000000000000000000000000000000000000..e3249f5f7f9253abb783d200c3469167ffa8d7ad
--- /dev/null
+++ b/lbmpy/methods/cumulantbased/fourth_order_correction.py
@@ -0,0 +1,237 @@
+import sympy as sp
+
+from lbmpy.moment_transforms import PRE_COLLISION_MONOMIAL_CUMULANT, POST_COLLISION_CUMULANT
+from lbmpy.moments import MOMENT_SYMBOLS, statistical_quantity_symbol
+from lbmpy.stencils import Stencil, LBStencil
+from pystencils import Assignment
+from pystencils.simp.assignment_collection import AssignmentCollection
+from .cumulantbasedmethod import CumulantBasedLbMethod
+
+X, Y, Z = MOMENT_SYMBOLS
+CORRECTED_FOURTH_ORDER_POLYNOMIALS = [X ** 2 * Y ** 2 - 2 * X ** 2 * Z ** 2 + Y ** 2 * Z ** 2,
+ X ** 2 * Y ** 2 + X ** 2 * Z ** 2 - 2 * Y ** 2 * Z ** 2,
+ X ** 2 * Y ** 2 + X ** 2 * Z ** 2 + Y ** 2 * Z ** 2,
+ X ** 2 * Y * Z,
+ X * Y ** 2 * Z,
+ X * Y * Z ** 2]
+FOURTH_ORDER_CORRECTION_SYMBOLS = sp.symbols("corr_fourth_:6")
+FOURTH_ORDER_RELAXATION_RATE_SYMBOLS = sp.symbols("corr_rr_:10")
+
+
+def add_fourth_order_correction(collision_rule, shear_relaxation_rate, bulk_relaxation_rate, limiter):
+ """Adds the fourth order correction terms (:cite:`geier2017`, eq. 44-49) to a given polynomial D3Q27
+ cumulant collision rule."""
+ method = collision_rule.method
+
+ if not isinstance(method, CumulantBasedLbMethod) or method.stencil != LBStencil(Stencil.D3Q27):
+ raise ValueError("Fourth-order correction is only defined for D3Q27 cumulant methods.")
+
+ polynomials = method.cumulants
+ rho = method.zeroth_order_equilibrium_moment_symbol
+
+ if not (set(CORRECTED_FOURTH_ORDER_POLYNOMIALS) < set(polynomials)):
+ raise ValueError("Fourth order correction requires polynomial cumulants "
+ f"{', '.join(CORRECTED_FOURTH_ORDER_POLYNOMIALS)} to be present")
+
+ # Call
+ # PC1 = (x^2 * y^2 - 2 * x^2 * z^2 + y^2 * z^2),
+ # PC2 = (x^2 * y^2 + x^2 * z^2 - 2 * y^2 * z^2),
+ # PC3 = (x^2 * y^2 + x^2 * z^2 + y^2 * z^2),
+ # PC4 = (x^2 * y * z),
+ # PC5 = (x * y^2 * z),
+ # PC6 = (x * y * z^2)
+ poly_symbols = [sp.Symbol(f'{POST_COLLISION_CUMULANT}_{polynomials.index(poly)}')
+ for poly in CORRECTED_FOURTH_ORDER_POLYNOMIALS]
+
+ a_symb, b_symb = sp.symbols("a_corr, b_corr")
+ a, b = get_optimal_additional_parameters(shear_relaxation_rate, bulk_relaxation_rate)
+ correction_terms = get_fourth_order_correction_terms(method.relaxation_rate_dict, rho, a_symb, b_symb)
+ optimal_parametrisation = get_optimal_parametrisation_with_limiters(collision_rule, shear_relaxation_rate,
+ bulk_relaxation_rate, limiter)
+
+ subexp_dict = collision_rule.subexpressions_dict
+ subexp_dict = {**subexp_dict,
+ a_symb: a,
+ b_symb: b,
+ **correction_terms.subexpressions_dict,
+ **optimal_parametrisation.subexpressions_dict,
+ **correction_terms.main_assignments_dict,
+ **optimal_parametrisation.main_assignments_dict}
+ for sym, cor in zip(poly_symbols, FOURTH_ORDER_CORRECTION_SYMBOLS):
+ subexp_dict[sym] += cor
+
+ collision_rule.set_sub_expressions_from_dict(subexp_dict)
+ collision_rule.topological_sort()
+
+ return collision_rule
+
+
+def get_optimal_additional_parameters(shear_relaxation_rate, bulk_relaxation_rate):
+ """
+ Calculates the optimal parametrization for the additional parameters provided in :cite:`geier2017`
+ Equations 115-116.
+ """
+
+ omega_1 = shear_relaxation_rate
+ omega_2 = bulk_relaxation_rate
+
+ omega_11 = omega_1 * omega_1
+ omega_12 = omega_1 * omega_2
+ omega_22 = omega_2 * omega_2
+
+ two = sp.Float(2)
+ three = sp.Float(3)
+ four = sp.Float(4)
+ nine = sp.Float(9)
+
+ denom_ab = (omega_1 - omega_2) * (omega_2 * (two + three * omega_1) - sp.Float(8) * omega_1)
+
+ a = (four * omega_11 + two * omega_12 * (omega_1 - sp.Float(6)) + omega_22 * (
+ omega_1 * (sp.Float(10) - three * omega_1) - four)) / denom_ab
+ b = (four * omega_12 * (nine * omega_1 - sp.Float(16)) - four * omega_11 - two * omega_22 * (
+ two + nine * omega_1 * (omega_1 - two))) / (three * denom_ab)
+
+ return a, b
+
+
+def get_fourth_order_correction_terms(rrate_dict, rho, a, b):
+ pc1, pc2, pc3, pc4, pc5, pc6 = CORRECTED_FOURTH_ORDER_POLYNOMIALS
+
+ omega_1 = rrate_dict[X * Y]
+ omega_2 = rrate_dict[X ** 2 + Y ** 2 + Z ** 2]
+
+ try:
+ omega_6 = rrate_dict[pc1]
+ assert omega_6 == rrate_dict[pc2], \
+ "Cumulants (x^2 - y^2) and (x^2 - z^2) must have the same relaxation rate"
+ omega_7 = rrate_dict[pc3]
+ omega_8 = rrate_dict[pc4]
+ assert omega_8 == rrate_dict[pc5] == rrate_dict[pc6]
+ except IndexError:
+ raise ValueError("For the fourth order correction, all six polynomial cumulants"
+ "(x^2 * y^2 - 2 * x^2 * z^2 + y^2 * z^2),"
+ "(x^2 * y^2 + x^2 * z^2 - 2 * y^2 * z^2),"
+ "(x^2 * y^2 + x^2 * z^2 + y^2 * z^2),"
+ "(x^2 * y * z), (x * y^2 * z) and (x * y * z^2) must be present!")
+
+ dxu, dyv, dzw, dxvdyu, dxwdzu, dywdzv = sp.symbols('Dx, Dy, Dz, DxvDyu, DxwDzu, DywDzv')
+ c_xy = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 1, 0))
+ c_xz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 0, 1))
+ c_yz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 1, 1))
+ c_xx = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (2, 0, 0))
+ c_yy = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 2, 0))
+ c_zz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 0, 2))
+
+ cor1, cor2, cor3, cor4, cor5, cor6 = FOURTH_ORDER_CORRECTION_SYMBOLS
+
+ one = sp.Float(1)
+ two = sp.Float(2)
+ three = sp.Float(3)
+
+ # Derivative Approximations
+ subexpressions = [
+ Assignment(dxu, - omega_1 / (two * rho) * (two * c_xx - c_yy - c_zz)
+ - omega_2 / (two * rho) * (c_xx + c_yy + c_zz - rho)),
+ Assignment(dyv, dxu + (three * omega_1) / (two * rho) * (c_xx - c_yy)),
+ Assignment(dzw, dxu + (three * omega_1) / (two * rho) * (c_xx - c_zz)),
+ Assignment(dxvdyu, - three * omega_1 / rho * c_xy),
+ Assignment(dxwdzu, - three * omega_1 / rho * c_xz),
+ Assignment(dywdzv, - three * omega_1 / rho * c_yz)]
+
+ one_half = sp.Rational(1, 2)
+ one_over_three = sp.Rational(1, 3)
+ two_over_three = sp.Rational(2, 3)
+ four_over_three = sp.Rational(4, 3)
+
+ # Correction Terms
+ main_assignments = [
+ Assignment(cor1, two_over_three * (one / omega_1 - one_half) * omega_6 * a * rho * (dxu - two * dyv + dzw)),
+ Assignment(cor2, two_over_three * (one / omega_1 - one_half) * omega_6 * a * rho * (dxu + dyv - two * dzw)),
+ Assignment(cor3, - four_over_three * (one / omega_1 - one_half) * omega_7 * a * rho * (dxu + dyv + dzw)),
+ Assignment(cor4, - one_over_three * (one / omega_1 - one_half) * omega_8 * b * rho * dywdzv),
+ Assignment(cor5, - one_over_three * (one / omega_1 - one_half) * omega_8 * b * rho * dxwdzu),
+ Assignment(cor6, - one_over_three * (one / omega_1 - one_half) * omega_8 * b * rho * dxvdyu)]
+
+ return AssignmentCollection(main_assignments=main_assignments, subexpressions=subexpressions)
+
+
+def get_optimal_parametrisation_with_limiters(collision_rule, shear_relaxation_rate, bulk_relaxation_rate, limiter):
+ """
+ Calculates the optimal parametrization for the relaxation rates provided in :cite:`geier2017`
+ Equations 112-114.
+ """
+
+ # if omega numbers
+ # assert omega_1 != omega_2, "Relaxation rates associated with shear and bulk viscosity must not be identical."
+ # assert omega_1 > 7/4
+
+ omega_1 = FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[0]
+ omega_2 = FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[1]
+
+ non_limited_omegas = sp.symbols("non_limited_omega_3:6")
+ limited_omegas = sp.symbols("limited_omega_3:6")
+
+ omega_11 = omega_1 * omega_1
+ omega_12 = omega_1 * omega_2
+ omega_22 = omega_2 * omega_2
+
+ one = sp.Float(1)
+ two = sp.Float(2)
+ three = sp.Float(3)
+ five = sp.Float(5)
+ six = sp.Float(6)
+ seven = sp.Float(7)
+ eight = sp.Float(8)
+ nine = sp.Float(9)
+
+ omega_3 = (eight * (omega_1 - two) * (omega_2 * (three * omega_1 - one) - five * omega_1)) / (
+ eight * (five - two * omega_1) * omega_1 + omega_2 * (eight + omega_1 * (nine * omega_1 - sp.Float(26))))
+
+ omega_4 = (eight * (omega_1 - two) * (omega_1 + omega_2 * (three * omega_1 - seven))) / (
+ omega_2 * (sp.Float(56) - sp.Float(42) * omega_1 + nine * omega_11) - eight * omega_1)
+
+ omega_5 = (sp.Float(24) * (omega_1 - two) * (sp.Float(4) * omega_11 + omega_12 * (
+ sp.Float(18) - sp.Float(13) * omega_1) + omega_22 * (two + omega_1 * (
+ six * omega_1 - sp.Float(11))))) / (sp.Float(16) * omega_11 * (omega_1 - six) - two * omega_12 * (
+ sp.Float(216) + five * omega_1 * (nine * omega_1 - sp.Float(46))) + omega_22 * (omega_1 * (
+ three * omega_1 - sp.Float(10)) * (sp.Float(15) * omega_1 - sp.Float(28)) - sp.Float(48)))
+
+ rho = collision_rule.method.zeroth_order_equilibrium_moment_symbol
+
+ # add limiters to improve stability
+ c_xyy = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 2, 0))
+ c_xzz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 0, 2))
+ c_xyz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 1, 1))
+ abs_c_xyy_plus_xzz = sp.Abs(c_xyy + c_xzz)
+ abs_c_xyy_minus_xzz = sp.Abs(c_xyy - c_xzz)
+ abs_c_xyz = sp.Abs(c_xyz)
+
+ limited_omega_3 = non_limited_omegas[0] + ((one - non_limited_omegas[0]) * abs_c_xyy_plus_xzz) / \
+ (rho * limiter + abs_c_xyy_plus_xzz)
+ limited_omega_4 = non_limited_omegas[1] + ((one - non_limited_omegas[1]) * abs_c_xyy_minus_xzz) / \
+ (rho * limiter + abs_c_xyy_minus_xzz)
+ limited_omega_5 = non_limited_omegas[2] + ((one - non_limited_omegas[2]) * abs_c_xyz) / (rho * limiter + abs_c_xyz)
+
+ subexpressions = [
+ Assignment(non_limited_omegas[0], omega_3),
+ Assignment(non_limited_omegas[1], omega_4),
+ Assignment(non_limited_omegas[2], omega_5),
+ Assignment(limited_omegas[0], limited_omega_3),
+ Assignment(limited_omegas[1], limited_omega_4),
+ Assignment(limited_omegas[2], limited_omega_5)]
+
+ # Correction Terms
+ main_assignments = [
+ Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[0], shear_relaxation_rate),
+ Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[1], bulk_relaxation_rate),
+ Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[2], limited_omegas[0]),
+ Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[3], limited_omegas[1]),
+ Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[4], limited_omegas[2]),
+ Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[5], one),
+ Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[6], one),
+ Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[7], one),
+ Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[8], one),
+ Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[9], one),
+ ]
+
+ return AssignmentCollection(main_assignments=main_assignments, subexpressions=subexpressions)
diff --git a/lbmpy/methods/momentbased/entropic.py b/lbmpy/methods/momentbased/entropic.py
index a9be8b1de0c07abb57b58c8074e7d7616a66c5ca..fae51497f82412cc7aa7468174d0e3eb424a7645 100644
--- a/lbmpy/methods/momentbased/entropic.py
+++ b/lbmpy/methods/momentbased/entropic.py
@@ -286,7 +286,7 @@ def _get_relaxation_rates(collision_rule):
omega_s = get_shear_relaxation_rate(method)
# if the shear relaxation rate is not specified as a symbol look for its symbolic counter part in the subs dict
- for symbolic_rr, rr in method.subs_dict_relxation_rate.items():
+ for symbolic_rr, rr in method.subs_dict_relaxation_rate.items():
if omega_s == rr:
omega_s = symbolic_rr
diff --git a/lbmpy/methods/momentbased/momentbasedmethod.py b/lbmpy/methods/momentbased/momentbasedmethod.py
index 2ce9d51580b6c9d8ea37344fe2d5e606feaae87f..79aee5c813b71c358e93964b52b84b6b1532269f 100644
--- a/lbmpy/methods/momentbased/momentbasedmethod.py
+++ b/lbmpy/methods/momentbased/momentbasedmethod.py
@@ -16,8 +16,8 @@ from lbmpy.moment_transforms import PdfsToMomentsByChimeraTransform
class MomentBasedLbMethod(AbstractLbMethod):
"""
Moment based LBM is a class to represent the single (SRT), two (TRT) and multi relaxation time (MRT) methods.
- These methods work by transforming the pdfs into moment space using a linear transformation. In the moment
- space each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These
+ These methods work by transforming the pdfs into moment space using a linear transformation. In the moment-space
+ each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These
equilibrium moments can e.g. be determined by taking the equilibrium moments of the continuous Maxwellian.
Parameters:
@@ -378,7 +378,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
subexpressions += m_post_to_f_post_eqs.subexpressions
main_assignments = m_post_to_f_post_eqs.main_assignments
else:
- # For SRT, TRT by default, and whenever customly required, derive equations entirely in
+ # For SRT, TRT by default, and whenever customarily required, derive equations entirely in
# population space
if self._zero_centered and not self._equilibrium.deviation_only:
raise Exception("Can only derive population-space equations for zero-centered storage"
diff --git a/lbmpy/scenarios.py b/lbmpy/scenarios.py
index 861bad743f5022af33983a132548b879c9bc5729..484b3ad6833eae588bcdd502ce7dfadf603e597b 100644
--- a/lbmpy/scenarios.py
+++ b/lbmpy/scenarios.py
@@ -42,10 +42,10 @@ def create_fully_periodic_flow(initial_velocity, periodicity_in_kernel=False, lb
initial_velocity: numpy array that defines an initial velocity for each cell. The shape of this
array determines the domain size.
periodicity_in_kernel: don't use boundary handling for periodicity, but directly generate the kernel periodic
- lbm_kernel: a LBM function, which would otherwise automatically created
+ lbm_kernel: an LBM function, which would otherwise automatically created
data_handling: data handling instance that is used to create the necessary fields. If a data handling is
passed the periodicity and parallel arguments are ignored.
- parallel: True for distributed memory parallelization with walberla
+ parallel: True for distributed memory parallelization with waLBerla
kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
Returns:
@@ -76,9 +76,9 @@ def create_lid_driven_cavity(domain_size=None, lid_velocity=0.005, lbm_kernel=No
Args:
domain_size: tuple specifying the number of cells in each dimension
lid_velocity: x velocity of lid in lattice coordinates.
- lbm_kernel: a LBM function, which would otherwise automatically created
+ lbm_kernel: an LBM function, which would otherwise automatically created
kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
- parallel: True for distributed memory parallelization with walberla
+ parallel: True for distributed memory parallelization with waLBerla
data_handling: see documentation of :func:`create_fully_periodic_flow`
Returns:
instance of :class:`Scenario`
diff --git a/lbmpy_tests/cumulantmethod/test_flow_around_sphere.py b/lbmpy_tests/cumulantmethod/test_flow_around_sphere.py
index 6db3fa9fbd3c91127b305cfda4864464763b2d64..ad5707bd13abcf14905455987aa46ce00dda4356 100644
--- a/lbmpy_tests/cumulantmethod/test_flow_around_sphere.py
+++ b/lbmpy_tests/cumulantmethod/test_flow_around_sphere.py
@@ -16,7 +16,7 @@ from pystencils import create_kernel, create_data_handling, Assignment, Target,
from pystencils.slicing import slice_from_direction, get_slice_before_ghost_layer
-def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
+def flow_around_sphere(stencil, galilean_correction, fourth_order_correction, L_LU, total_steps):
if galilean_correction and stencil.Q != 27:
pytest.skip()
@@ -38,7 +38,7 @@ def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
lbm_config = LBMConfig(stencil=stencil, method=Method.CUMULANT, relaxation_rate=omega_v,
compressible=True,
- galilean_correction=galilean_correction)
+ galilean_correction=galilean_correction, fourth_order_correction=fourth_order_correction)
lbm_opt = LBMOptimisation(pre_simplification=True)
config = CreateKernelConfig(target=target)
@@ -135,14 +135,22 @@ def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
@pytest.mark.parametrize('galilean_correction', [False, True])
-def test_flow_around_sphere_short(stencil, galilean_correction):
+@pytest.mark.parametrize('fourth_order_correction', [False, True])
+def test_flow_around_sphere_short(stencil, galilean_correction, fourth_order_correction):
pytest.importorskip('cupy')
- flow_around_sphere(LBStencil(stencil), galilean_correction, 5, 200)
+ stencil = LBStencil(stencil)
+ if fourth_order_correction and stencil.Q != 27:
+ pytest.skip("Fourth-order correction only defined for D3Q27 stencil.")
+ flow_around_sphere(stencil, galilean_correction, fourth_order_correction, 5, 200)
@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
@pytest.mark.parametrize('galilean_correction', [False, True])
+@pytest.mark.parametrize('fourth_order_correction', [False, 0.01])
@pytest.mark.longrun
-def test_flow_around_sphere_long(stencil, galilean_correction):
+def test_flow_around_sphere_long(stencil, galilean_correction, fourth_order_correction):
pytest.importorskip('cupy')
- flow_around_sphere(LBStencil(stencil), galilean_correction, 20, 3000)
+ stencil = LBStencil(stencil)
+ if fourth_order_correction and stencil.Q != 27:
+ pytest.skip("Fourth-order correction only defined for D3Q27 stencil.")
+ flow_around_sphere(stencil, galilean_correction, fourth_order_correction, 20, 3000)
diff --git a/lbmpy_tests/full_scenarios/shear_wave/scenario_shear_wave.py b/lbmpy_tests/full_scenarios/shear_wave/scenario_shear_wave.py
index cd5e73b7953216c3dd22c79490ecf29c178dc944..8d91315eb447510f336082a64448435eab068731 100644
--- a/lbmpy_tests/full_scenarios/shear_wave/scenario_shear_wave.py
+++ b/lbmpy_tests/full_scenarios/shear_wave/scenario_shear_wave.py
@@ -5,15 +5,16 @@
Chapter 5.1
"""
import numpy as np
-import sympy as sp
-
import pytest
-from pystencils import Target
+import sympy as sp
-from lbmpy.creationfunctions import update_with_default_parameters
+from lbmpy import LatticeBoltzmannStep, LBStencil
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation
+from lbmpy.db import LbmpyJsonSerializer
+from lbmpy.enums import Method, Stencil
from lbmpy.relaxationrates import (
relaxation_rate_from_lattice_viscosity, relaxation_rate_from_magic_number)
-from lbmpy.scenarios import create_fully_periodic_flow
+from pystencils import Target, create_data_handling, CreateKernelConfig
def get_exponent_term(l, **kwargs):
@@ -53,14 +54,13 @@ def get_analytical_solution(l, l_0, u_0, v_0, nu, t):
def plot_y_velocity(vel, **kwargs):
import matplotlib.pyplot as plt
- from matplotlib import cm
- vel = vel[:, vel.shape[1]//2, :, 1]
+ vel = vel[:, vel.shape[1] // 2, :, 1]
ranges = [np.arange(n, dtype=float) for n in vel.shape]
x, y = np.meshgrid(*ranges, indexing='ij')
fig = plt.gcf()
ax = fig.gca(projection='3d')
- ax.plot_surface(x, y, vel, cmap=cm.coolwarm, rstride=2, cstride=2,
+ ax.plot_surface(x, y, vel, cmap='coolwarm', rstride=2, cstride=2,
linewidth=0, antialiased=True, **kwargs)
@@ -71,9 +71,9 @@ def fit_and_get_slope(x_values, y_values):
def get_phase_error(phases, evaluation_interval):
- xValues = np.arange(len(phases) * evaluation_interval, step=evaluation_interval)
- phaseError = fit_and_get_slope(xValues, phases)
- return phaseError
+ x_values = np.arange(len(phases) * evaluation_interval, step=evaluation_interval)
+ phase_error = fit_and_get_slope(x_values, phases)
+ return phase_error
def get_viscosity(abs_values, evaluation_interval, l):
@@ -92,33 +92,44 @@ def get_amplitude_and_phase(vel_slice):
return fft_abs[max_index], fft_phase[max_index]
-def run(l, l_0, u_0, v_0, nu, y_size, periodicity_in_kernel, **kwargs):
+def run(l, l_0, u_0, v_0, nu, y_size, lbm_config, lbm_optimisation, config):
inv_result = {
'viscosity_error': np.nan,
'phase_error': np.nan,
'mlups': np.nan,
}
- if 'initial_velocity' in kwargs:
- del kwargs['initial_velocity']
+ if lbm_config.initial_velocity:
+ # no manually specified initial velocity allowed in shear wave setup
+ lbm_config.initial_velocity = None
- print("Running size l=%d nu=%f " % (l, nu), kwargs)
+ print(f"Running size l={l} nu={nu}")
initial_vel_arr = get_initial_velocity_field(l, l_0, u_0, v_0, y_size)
omega = relaxation_rate_from_lattice_viscosity(nu)
- kwargs['relaxation_rates'] = [sp.sympify(r) for r in kwargs['relaxation_rates']]
- if 'entropic' not in kwargs or not kwargs['entropic']:
- kwargs['relaxation_rates'] = [r.subs(sp.Symbol("omega"), omega) for r in kwargs['relaxation_rates']]
+ if lbm_config.fourth_order_correction and omega < 1.75:
+ pytest.skip("Fourth-order correction only allowed for omega >= 1.75.")
+
+ lbm_config.relaxation_rates = [sp.sympify(r) for r in lbm_config.relaxation_rates]
+ lbm_config.relaxation_rates = [r.subs(sp.Symbol("omega"), omega) for r in lbm_config.relaxation_rates]
+
+ periodicity_in_kernel = (lbm_optimisation.builtin_periodicity != (False, False, False))
+ domain_size = initial_vel_arr.shape[:-1]
- scenario = create_fully_periodic_flow(initial_vel_arr, periodicity_in_kernel=periodicity_in_kernel, **kwargs)
+ data_handling = create_data_handling(domain_size, periodicity=not periodicity_in_kernel,
+ default_ghost_layers=1, parallel=False)
- if 'entropic' in kwargs and kwargs['entropic']:
- scenario.kernelParams['omega'] = kwargs['relaxation_rates'][0].subs(sp.Symbol("omega"), omega)
+ scenario = LatticeBoltzmannStep(data_handling=data_handling, name="periodic_scenario", lbm_config=lbm_config,
+ lbm_optimisation=lbm_optimisation, config=config)
+ for b in scenario.data_handling.iterate(ghost_layers=False):
+ np.copyto(b[scenario.velocity_data_name], initial_vel_arr[b.global_slice])
+ scenario.set_pdf_fields_from_macroscopic_values()
total_time_steps = 20000 * (l // l_0) ** 2
initial_time_steps = 11000 * (l // l_0) ** 2
eval_interval = 1000 * (l // l_0) ** 2
scenario.run(initial_time_steps)
if np.isnan(scenario.velocity_slice()).any():
+ print(" Result", inv_result)
return inv_result
magnitudes = []
@@ -149,65 +160,77 @@ def run(l, l_0, u_0, v_0, nu, y_size, periodicity_in_kernel, **kwargs):
def create_full_parameter_study():
from pystencils.runhelper import ParameterStudy
- params = {
+ setup_params = {
'l_0': 32,
'u_0': 0.096,
'v_0': 0.1,
- 'ySize': 1,
- 'periodicityInKernel': True,
- 'stencil': 'D3Q27',
- 'compressible': True,
+ 'y_size': 1
}
+
+ omega, omega_f = sp.symbols("omega, omega_f")
+
ls = [32 * 2 ** i for i in range(0, 5)]
nus = [1e-2, 1e-3, 1e-4, 1e-5]
- srt_and_trt_methods = [{'method': method,
- 'stencil': stencil,
- 'compressible': comp,
- 'relaxation_rates': ["omega", str(relaxation_rate_from_magic_number(sp.Symbol("omega")))],
- 'equilibrium_order': eqOrder,
- 'continuous_equilibrium': mbEq,
- 'optimization': {'target': Target.CPU, 'split': True, 'cse_pdfs': True}}
- for method in ('srt', 'trt')
- for stencil in ('D3Q19', 'D3Q27')
+ srt_and_trt_methods = [LBMConfig(method=method,
+ stencil=LBStencil(stencil),
+ compressible=comp,
+ relaxation_rates=[omega, str(relaxation_rate_from_magic_number(omega))],
+ equilibrium_order=eqOrder,
+ continuous_equilibrium=mbEq)
+ for method in (Method.SRT, Method.TRT)
+ for stencil in (Stencil.D3Q19, Stencil.D3Q27)
for comp in (True, False)
for eqOrder in (1, 2, 3)
for mbEq in (True, False)]
- methods = [{'method': 'trt', 'relaxation_rates': ["omega", 1]},
- {'method': 'mrt3', 'relaxation_rates': ["omega", 1, 1], 'equilibrium_order': 2},
- {'method': 'mrt3', 'relaxation_rates': ["omega", 1, 1], 'equilibrium_order': 3},
- {'method': 'mrt3', 'cumulant': True, 'relaxation_rates': ["omega", 1, 1], # cumulant
- 'continuous_equilibrium': True, 'equilibrium_order': 3,
- 'optimization': {'cse_global': True}},
- {'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True, # entropic order 2
- 'continuous_equilibrium': True, 'equilibrium_order': 2},
- {'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True, # entropic order 3
- 'continuous_equilibrium': True, 'equilibrium_order': 3},
-
- # entropic cumulant:
- {'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,
- 'cumulant': True, 'continuous_equilibrium': True, 'equilibrium_order': 3},
- {'method': 'trt-kbc-n2', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,
- 'cumulant': True, 'continuous_equilibrium': True, 'equilibrium_order': 3},
- {'method': 'mrt3', 'relaxation_rates': ["omega", "omega_f", "omega_f"], 'entropic': True,
- 'cumulant': True, 'continuous_equilibrium': True, 'equilibrium_order': 3},
+ optimization_srt_trt = LBMOptimisation(split=True, cse_pdfs=True)
+
+ config = CreateKernelConfig(target=Target.CPU)
+
+ methods = [LBMConfig(method=Method.TRT, relaxation_rates=[omega, 1]),
+ LBMConfig(method=Method.MRT, relaxation_rates=[omega], equilibrium_order=2),
+ LBMConfig(method=Method.MRT, relaxation_rates=[omega], equilibrium_order=3),
+ LBMConfig(method=Method.CUMULANT, relaxation_rates=[omega], # cumulant
+ compressible=True, continuous_equilibrium=True, equilibrium_order=3),
+ LBMConfig(method=Method.CUMULANT, relaxation_rates=[omega], # cumulant with fourth-order correction
+ compressible=True, continuous_equilibrium=True, fourth_order_correction=0.1),
+ LBMConfig(method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f], entropic=True,
+ zero_centered=False, # entropic order 2
+ continuous_equilibrium=True, equilibrium_order=2),
+ LBMConfig(method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f], entropic=True,
+ zero_centered=False, # entropic order 3
+ continuous_equilibrium=True, equilibrium_order=3),
+
+ # entropic cumulant: not supported for the moment
+ # LBMConfig(method=Method.CUMULANT, relaxation_rates=["omega", "omega_f"], entropic=True, zero_centered=False,
+ # compressible=True, continuous_equilibrium=True, equilibrium_order=3)
]
- parameter_study = ParameterStudy(run, database_connector="shear_wave_db")
+ parameter_study = ParameterStudy(run, database_connector="shear_wave_db",
+ serializer_info=('lbmpy_serializer', LbmpyJsonSerializer))
for L in ls:
for nu in nus:
- for methodParams in methods + srt_and_trt_methods:
- simulation_params = params.copy()
- simulation_params.update(methodParams)
- if 'optimization' not in simulation_params:
- simulation_params['optimization'] = {}
- new_params, new_opt_params = update_with_default_parameters(simulation_params,
- simulation_params['optimization'], False)
- simulation_params = new_params
- simulation_params['optimization'] = new_opt_params
-
- simulation_params['L'] = L
+ for methodParams in srt_and_trt_methods:
+ simulation_params = setup_params.copy()
+
+ simulation_params['lbm_config'] = methodParams
+ simulation_params['lbm_optimisation'] = optimization_srt_trt
+ simulation_params['config'] = config
+
+ simulation_params['l'] = L
+ simulation_params['nu'] = nu
+ l_0 = simulation_params['l_0']
+ parameter_study.add_run(simulation_params.copy(), weight=(L / l_0) ** 4)
+
+ for methodParams in methods:
+ simulation_params = setup_params.copy()
+
+ simulation_params['lbm_config'] = methodParams
+ simulation_params['lbm_optimisation'] = LBMOptimisation()
+ simulation_params['config'] = config
+
+ simulation_params['l'] = L
simulation_params['nu'] = nu
l_0 = simulation_params['l_0']
parameter_study.add_run(simulation_params.copy(), weight=(L / l_0) ** 4)
@@ -217,13 +240,23 @@ def create_full_parameter_study():
def test_shear_wave():
pytest.importorskip('cupy')
params = {
+ 'y_size': 1,
'l_0': 32,
'u_0': 0.096,
'v_0': 0.1,
- 'stencil': 'D3Q19',
- 'compressible': True,
- "optimization": {"target": Target.GPU}
+ 'l': 32,
+ 'nu': 1e-2,
}
- run(32, nu=1e-2, equilibrium_order=2, method='srt', y_size=1, periodicity_in_kernel=True,
- relaxation_rates=[sp.Symbol("omega"), 5, 5], continuous_equilibrium=True, **params)
+
+ kernel_config = CreateKernelConfig(target=Target.GPU)
+ lbm_config = LBMConfig(method=Method.SRT, relaxation_rates=[sp.Symbol("omega")], stencil=LBStencil(Stencil.D3Q27),
+ compressible=True, continuous_equilibrium=True, equilibrium_order=2)
+
+ run(**params, lbm_config=lbm_config, config=kernel_config, lbm_optimisation=LBMOptimisation())
+
+
+@pytest.mark.longrun
+def test_all_scenarios():
+ parameter_study = create_full_parameter_study()
+ parameter_study.run()
diff --git a/lbmpy_tests/test_force.py b/lbmpy_tests/test_force.py
index 7ac94f42b5cf649147cbefa5830d2221c628cc14..f0f59bb5490cb753d4f7b0e5efb70ec65c2b9caf 100644
--- a/lbmpy_tests/test_force.py
+++ b/lbmpy_tests/test_force.py
@@ -197,7 +197,7 @@ def test_literature(force_model, stencil, method):
assert len(omega_momentum) == 1
omega_momentum = omega_momentum[0]
- subs_dict = lb_method.subs_dict_relxation_rate
+ subs_dict = lb_method.subs_dict_relaxation_rate
force_term = sp.simplify(lb_method.force_model(lb_method).subs(subs_dict))
u = sp.Matrix(lb_method.first_order_equilibrium_moment_symbols)
rho = lb_method.conserved_quantity_computation.density_symbol
@@ -282,7 +282,7 @@ def test_modes_central_moment(force_model, compressible):
compressible=compressible, force_model=force_model, force=tuple(F))
method = create_lb_method(lbm_config=lbm_config)
- subs_dict = method.subs_dict_relxation_rate
+ subs_dict = method.subs_dict_relaxation_rate
force_moments = method.force_model.central_moment_space_forcing(method)
force_moments = force_moments.subs(subs_dict)
@@ -310,7 +310,7 @@ def test_symmetric_forcing_equivalence(force_model, compressible):
if not method.force_model.has_symmetric_central_moment_forcing:
return True
- subs_dict = method.subs_dict_relxation_rate
+ subs_dict = method.subs_dict_relaxation_rate
force_moments = method.force_model.central_moment_space_forcing(method)
force_moments = force_moments.subs(subs_dict)
@@ -335,7 +335,7 @@ def test_modes_central_moment_longrun(stencil, force_model, compressible):
compressible=compressible, force_model=force_model, force=tuple(F))
method = create_lb_method(lbm_config=lbm_config)
- subs_dict = method.subs_dict_relxation_rate
+ subs_dict = method.subs_dict_relaxation_rate
force_moments = method.force_model.moment_space_forcing(method)
force_moments = force_moments.subs(subs_dict)
@@ -359,7 +359,7 @@ def _check_modes(stencil, force_model, compressible):
compressible=compressible, force_model=force_model, force=tuple(F))
method = create_lb_method(lbm_config=lbm_config)
- subs_dict = method.subs_dict_relxation_rate
+ subs_dict = method.subs_dict_relaxation_rate
force_moments = method.force_model.moment_space_forcing(method)
force_moments = force_moments.subs(subs_dict)
diff --git a/lbmpy_tests/test_json_serializer.py b/lbmpy_tests/test_json_serializer.py
new file mode 100644
index 0000000000000000000000000000000000000000..3429f12bb6ba0eff62e76e127707f08b290d8547
--- /dev/null
+++ b/lbmpy_tests/test_json_serializer.py
@@ -0,0 +1,40 @@
+"""
+Test the lbmpy-specific JSON encoder and serializer as used in the Database class.
+"""
+
+import tempfile
+import pystencils as ps
+
+from lbmpy import Stencil, Method, ForceModel
+from lbmpy.advanced_streaming import Timestep
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, LBStencil
+from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor
+
+from pystencils.runhelper import Database
+from lbmpy.db import LbmpyJsonSerializer
+
+
+def test_json_serializer():
+
+ stencil = LBStencil(Stencil.D3Q27)
+ q = stencil.Q
+ pdfs, pdfs_tmp = ps.fields(f"pdfs({q}), pdfs_tmp({q}): double[3D]", layout='fzyx')
+ density = ps.fields(f"rho: double[3D]", layout='fzyx')
+
+ from lbmpy.non_newtonian_models import CassonsParameters
+ cassons_params = CassonsParameters(0.2)
+
+ # create dummy lbmpy config
+ lbm_config = LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, force_model=ForceModel.GUO,
+ compressible=True, relaxation_rate=1.999, smagorinsky=True, galilean_correction=True,
+ cassons=cassons_params, density_input=density, kernel_type=StreamPullTwoFieldsAccessor,
+ timestep=Timestep.BOTH)
+
+ lbm_optimization = LBMOptimisation(cse_pdfs=False, cse_global=False, builtin_periodicity=(True, False, False),
+ symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp)
+
+ # create dummy database
+ temp_dir = tempfile.TemporaryDirectory()
+ db = Database(file=temp_dir.name, serializer_info=('lbmpy_serializer', LbmpyJsonSerializer))
+
+ db.save(params={'lbm_config': lbm_config, 'lbm_optimization': lbm_optimization}, result={'test': 'dummy'})
diff --git a/lbmpy_tests/test_shear_flow.py b/lbmpy_tests/test_shear_flow.py
index 9b22d0a5b0a8cd36ce53b0e32c649428959ecf9a..4f15d43ab5ae57fa5f3993b16be8e75fd5e357d6 100644
--- a/lbmpy_tests/test_shear_flow.py
+++ b/lbmpy_tests/test_shear_flow.py
@@ -1,5 +1,5 @@
"""
-Test shear flow velocity and pressureagainst analytical solutions
+Test shear flow velocity and pressure against analytical solutions
"""
@@ -205,7 +205,7 @@ def test_shear_flow(target, stencil_name, zero_centered):
p_expected = rho_0 * np.identity(dh.dim) / 3.0 + dynamic_viscosity * shear_rate / correction_factor * (
np.outer(shear_plane_normal, shear_direction) + np.transpose(np.outer(shear_plane_normal, shear_direction)))
- # Sustract the tensorproduct of the velosity to get the pressure
+ # Subtract the tensor product of the velocity to get the pressure
pressure_profile[:, 0, 0] -= rho_0 * profile[:, 0]**2
np.testing.assert_allclose(profile[:, 0], expected[1:-1], atol=1E-9)