Commit 21b1d786 authored by Helen Schottenhamml's avatar Helen Schottenhamml
Browse files

Merge branch 'NonNewtonian' into 'master'


See merge request pycodegen/lbmpy!118
parents c004c33a ae65dd25
Pipeline #39184 passed with stages
in 125 minutes and 17 seconds
......@@ -131,7 +131,7 @@ minimal-windows:
- env
- pip list
- python -c "import numpy"
- py.test -v -n $NUM_CORES -m "not (notebook or longrun)"
- py.test -v -m "not (notebook or longrun)"
stage: test
......@@ -792,7 +792,7 @@
"hash": "ca06c80c4febc35b85e85156d391051f9f4a8895eee3f708eb1f33a09d8697a0"
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
......@@ -806,7 +806,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.2"
"version": "3.9.9"
"nbformat": 4,
This diff is collapsed.
......@@ -239,4 +239,13 @@ journal = {Communications in Computational Physics}
publisher = {Elsevier {BV}},
author = {R. Ouared and B. Chopard},
journal = {Journal of Statistical Physics},
title = {Lattice {Boltzmann} Simulations of Blood Flow: {Non-Newtonian} Rheology and Clotting Processes},
year = {2005},
doi = {10.1007/s10955-005-8415-x},
publisher = {Springer Link},
@Comment{jabref-meta: databaseType:bibtex;}
......@@ -18,6 +18,7 @@ You can open the notebooks directly to play around with the code examples.
......@@ -63,6 +63,7 @@ import lbmpy.forcemodels as forcemodels
import lbmpy.methods.centeredcumulant.force_model as cumulant_force_model
from lbmpy.fieldaccess import CollideOnlyInplaceAccessor, PdfFieldAccessor, PeriodicTwoFieldsAccessor
from lbmpy.fluctuatinglb import add_fluctuations_to_collision_rule
from lbmpy.non_newtonian_models import add_cassons_model, CassonsParameters
from lbmpy.methods import (create_mrt_orthogonal, create_mrt_raw, create_central_moment,
create_srt, create_trt, create_trt_kbc)
from lbmpy.methods.creationfunctions import CollisionSpaceInfo
......@@ -225,6 +226,11 @@ class LBMConfig:
set to Smagorinsky constant to activate turbulence model, ``omega_output_field`` can be set to
write out adapted relaxation rates. If set to `True`, 0.12 is used as default smagorinsky constant.
cassons: CassonsParameters = False
Adds the Cassons model according to
The parameters are set with the ``CassonsParameters`` dataclass.
fluctuating: dict = False
Enables fluctuating lattice Boltzmann by randomizing collision process.
......@@ -628,8 +634,8 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
collision_rule = lb_method.get_collision_rule(pre_simplification=pre_simplification)
if lbm_config.entropic:
if lbm_config.smagorinsky:
raise ValueError("Choose either entropic or smagorinsky")
if lbm_config.smagorinsky or lbm_config.cassons:
raise ValueError("Choose either entropic, smagorinsky or cassons")
if lbm_config.entropic_newton_iterations:
if isinstance(lbm_config.entropic_newton_iterations, bool):
iterations = 3
......@@ -640,12 +646,18 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
collision_rule = add_entropy_condition(collision_rule, omega_output_field=lbm_config.omega_output_field)
elif lbm_config.smagorinsky:
if lbm_config.cassons:
raise ValueError("Cassons model can not be combined with Smagorinsky model")
smagorinsky_constant = 0.12 if lbm_config.smagorinsky is True else lbm_config.smagorinsky
collision_rule = add_smagorinsky_model(collision_rule, smagorinsky_constant,
if 'split_groups' in collision_rule.simplification_hints:
elif lbm_config.cassons:
collision_rule = add_cassons_model(collision_rule, parameter=lbm_config.cassons,
if lbm_config.output:
cqc = lb_method.conserved_quantity_computation
output_eqs = cqc.output_equations_from_pdfs(lb_method.pre_collision_pdf_symbols, lbm_config.output)
from dataclasses import dataclass
import sympy as sp
from lbmpy.relaxationrates import get_shear_relaxation_rate, lattice_viscosity_from_relaxation_rate, \
from lbmpy.utils import extract_shear_relaxation_rate, frobenius_norm, second_order_moment_tensor
from pystencils import Assignment
class CassonsParameters:
yield_stress: float
The yield stress controls the intensity of non-Newtonian behavior.
omega_min: float = 0.2
The minimal shear relaxation rate that is used as a lower bound
omega_max: float = 1.98
The maximal shear relaxation rate that is used as an upper bound
def add_cassons_model(collision_rule, parameter: CassonsParameters, omega_output_field=None):
r""" Adds the Cassons model to a lattice Boltzmann collision rule that can be found here :cite:`Casson`
The only parameter of the model is the so-called yield_stress. The main idea is that no strain rate is
observed below some stress. However, this leads to the problem that the modified relaxation rate might no longer
lead to stable LBM simulations. Thus, an upper and lower limit for the shear relaxation rate must be given.
All the parameters are combined in the `CassonsParameters` dataclass
yield_stress = parameter.yield_stress
omega_min = parameter.omega_min
omega_max = parameter.omega_max
method = collision_rule.method
equilibrium = method.equilibrium_distribution
omega_s = get_shear_relaxation_rate(method)
omega_s, found_symbolic_shear_relaxation = extract_shear_relaxation_rate(collision_rule, omega_s)
if not found_symbolic_shear_relaxation:
raise ValueError("For the Cassons model the shear relaxation rate has to be a symbol or it has to be "
"assigned to a single equation in the assignment list")
sigma, theta, rhs, mu, tau, adapted_omega = sp.symbols("sigma theta rhs mu tau omega_new")
rho = equilibrium.density if equilibrium.compressible else equilibrium.background_density
f_neq = sp.Matrix(method.pre_collision_pdf_symbols) - method.get_equilibrium_terms()
second_invariant_strain_rate_tensor = frobenius_norm(second_order_moment_tensor(f_neq, method.stencil))
eta = lattice_viscosity_from_relaxation_rate(omega_s)
one = sp.Rational(1, 1)
# rhs of equation 14 in
# Note that C_2 / C_4 = 3 for all configurations thus we directly insert it here
eq14 = one / (one - theta) * (one + sp.sqrt(theta * (one + rho / eta * sp.Rational(3, 2) * (one - theta))))
new_omega = one / tau
omega_cond = sp.Piecewise((omega_min, new_omega < omega_min), (omega_max, new_omega > omega_max), (new_omega, True))
eqs = [Assignment(sigma, second_invariant_strain_rate_tensor),
Assignment(theta, yield_stress / sigma),
Assignment(rhs, eq14),
Assignment(mu, eta * rhs ** 2),
Assignment(tau, one / relaxation_rate_from_lattice_viscosity(mu)),
Assignment(adapted_omega, omega_cond)]
collision_rule = collision_rule.new_with_substitutions({omega_s: adapted_omega}, substitute_on_lhs=False)
collision_rule.subexpressions += eqs
collision_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
if omega_output_field:
collision_rule.main_assignments.append(Assignment(, adapted_omega))
return collision_rule
import sympy as sp
from lbmpy.relaxationrates import get_shear_relaxation_rate
from lbmpy.utils import frobenius_norm, second_order_moment_tensor
from lbmpy.utils import extract_shear_relaxation_rate, frobenius_norm, second_order_moment_tensor
from pystencils import Assignment
......@@ -31,14 +31,7 @@ def add_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_fie
method = collision_rule.method
omega_s = get_shear_relaxation_rate(method)
found_symbolic_shear_relaxation = True
if isinstance(omega_s, float) or isinstance(omega_s, int):
found_symbolic_shear_relaxation = False
for eq in collision_rule.all_assignments:
if eq.rhs == omega_s:
found_symbolic_shear_relaxation = True
omega_s = eq.lhs
omega_s, found_symbolic_shear_relaxation = extract_shear_relaxation_rate(collision_rule, omega_s)
if not found_symbolic_shear_relaxation:
raise ValueError("For the smagorinsky model the shear relaxation rate has to be a symbol or it has to be "
......@@ -11,3 +11,18 @@ def frobenius_norm(matrix, factor=1):
"""Computes the Frobenius norm of a matrix defined as the square root of the sum of squared matrix elements
The optional factor is added inside the square root"""
return sp.sqrt(sum(i * i for i in matrix) * factor)
def extract_shear_relaxation_rate(collision_rule, shear_relaxation_rate):
"""Searches for the shear relaxation rate in the collision equations.
If the shear relaxation rate is assigned to a single assignment its lhs is returned.
Otherwise, the shear relaxation rate has to be a sympy symbol, or it can not be extracted"""
found_symbolic_shear_relaxation = True
if isinstance(shear_relaxation_rate, (float, int)):
found_symbolic_shear_relaxation = False
for eq in collision_rule.all_assignments:
if eq.rhs == shear_relaxation_rate:
found_symbolic_shear_relaxation = True
shear_relaxation_rate = eq.lhs
return shear_relaxation_rate, found_symbolic_shear_relaxation
......@@ -27,8 +27,8 @@ def test_creation(method_enum, double_precision):
@pytest.mark.parametrize('double_precision', [False, True])
@pytest.mark.parametrize('method_enum', [Method.SRT, Method.CENTRAL_MOMENT, Method.CUMULANT])
def test_scenario(method_enum, double_precision):
lbm_config = LBMConfig(method=method_enum, relaxation_rate=1.5, compressible=True)
config = ps.CreateKernelConfig(data_type="double" if double_precision else "float32")
lbm_config = LBMConfig(method=method_enum, relaxation_rate=1.45, compressible=True)
config = ps.CreateKernelConfig(data_type="float64" if double_precision else "float32")
sc = create_lid_driven_cavity((16, 16, 8), lbm_config=lbm_config, config=config)
code = ps.get_code_str(sc.ast)
......@@ -91,7 +91,7 @@ def test_advanced_initialization():
init_vel[:, height // 3: height // 3 * 2, 0] = -velocity_magnitude
# small random y velocity component
init_vel[:, :, 1] = 0.1 * velocity_magnitude * np.random.rand(width, height)
shear_flow_scenario = create_fully_periodic_flow(initial_velocity=init_vel, relaxation_rate=1.95)
shear_flow_scenario = create_fully_periodic_flow(initial_velocity=init_vel, relaxation_rate=1.99)
with pytest.raises(ValueError) as e:
shear_flow_scenario.run_iterative_initialization(max_steps=20000, check_residuum_after=500)
assert 'did not converge' in str(e.value)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment