From 223d1abd814edaac2bfa006ac30cb2d869cf2741 Mon Sep 17 00:00:00 2001
From: Rudolf Weeber <weeber@icp.uni-stuttgart.de>
Date: Thu, 20 Mar 2025 16:47:02 +0100
Subject: [PATCH] Fix fluctuating LB with zero centered pdf storage
---
src/lbmpy/fluctuatinglb.py | 2 +-
tests/test_fluctuating_lb.py | 40 ++++++++++++++++++++----------------
2 files changed, 23 insertions(+), 19 deletions(-)
diff --git a/src/lbmpy/fluctuatinglb.py b/src/lbmpy/fluctuatinglb.py
index cd035fc1..d0232775 100644
--- a/src/lbmpy/fluctuatinglb.py
+++ b/src/lbmpy/fluctuatinglb.py
@@ -42,7 +42,7 @@ def fluctuation_amplitude_from_temperature(method, temperature, c_s_sq=sp.Symbol
"""Produces amplitude equations according to (2.60) and (3.54) in Schiller08"""
normalization_factors = sp.matrix_multiply_elementwise(method.moment_matrix, method.moment_matrix) * \
sp.Matrix(method.weights)
- density = method.zeroth_order_equilibrium_moment_symbol
+ density = method._cqc.density_symbol
mu = temperature * density / c_s_sq
return [sp.sqrt(mu * norm * (1 - (1 - rr) ** 2))
for norm, rr in zip(normalization_factors, method.relaxation_rates)]
diff --git a/tests/test_fluctuating_lb.py b/tests/test_fluctuating_lb.py
index 7086898a..9fe76777 100644
--- a/tests/test_fluctuating_lb.py
+++ b/tests/test_fluctuating_lb.py
@@ -14,7 +14,7 @@ from pystencils.enums import Target
from pystencils.rng import PhiloxTwoDoubles
from lbmpy.creationfunctions import *
-from lbmpy.forcemodels import Guo
+from lbmpy.forcemodels import Schiller
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
import numpy as np
from lbmpy.enums import Stencil
@@ -137,7 +137,7 @@ def get_fluctuating_lb(
weighted=True,
zero_centered=zero_centered,
relaxation_rates=rr_getter,
- force_model=Guo(force=force_field.center_vector),
+ force_model=Schiller(force=force_field.center_vector),
fluctuating={"temperature": kT},
kernel_type="collide_only",
)
@@ -175,7 +175,7 @@ def get_fluctuating_lb(
collision.method,
velocity=(0,) * dh.dim,
pdfs=src.center_vector,
- density=rho.center,
+ density=rho_0
)
init_kernel = ps.create_kernel(init, ghost_layers=0).compile()
@@ -216,7 +216,7 @@ def get_fluctuating_lb(
"domain_size", [8, 60]
)
def test_resting_fluid(zero_centered: bool, domain_size: int, target=Target.CPU):
- rho_0 = 0.86
+ rho_0 = 0.86
kT = 4e-4
L = [domain_size] * 3
dh, time_loop = get_fluctuating_lb(
@@ -243,16 +243,18 @@ def test_resting_fluid(zero_centered: bool, domain_size: int, target=Target.CPU)
res_u = dh.gather_array("u").reshape((-1, 3))
res_rho = dh.gather_array("rho").reshape((-1,))
- # mass conservation
- np.testing.assert_allclose(np.mean(res_rho), rho_0, atol=3e-12)
+ # mass conservationo
+ # density per cell fluctuates, but toal mass is conserved
+ np.testing.assert_allclose(np.mean(res_rho), rho_0, atol=3E-12)
# momentum conservation
momentum = np.dot(res_rho, res_u)
np.testing.assert_allclose(momentum, [0, 0, 0], atol=1e-10)
- # temperature
+ # temperature (fluctuates around pre-set kT)
kinetic_energy = 1 / 2 * np.dot(res_rho, res_u * res_u) / np.prod(L)
- np.testing.assert_allclose(kinetic_energy, [kT / 2] * 3, atol=kT * 0.01)
+ kT_tol = 0.075 *(16/domain_size)**(3/2)
+ np.testing.assert_allclose(kinetic_energy, [kT / 2] * 3, rtol=kT_tol)
# velocity distribution
v_hist, v_bins = np.histogram(
@@ -271,12 +273,12 @@ def test_resting_fluid(zero_centered: bool, domain_size: int, target=Target.CPU)
v_expected.append(res)
v_expected = np.array(v_expected)
- # 10% accuracy on the entire histogram
- np.testing.assert_allclose(v_hist, v_expected, rtol=0.1)
- # 1% accuracy on the middle part
+ hist_tol_all = 0.75 *(16/domain_size)**(3/2)
+ np.testing.assert_allclose(v_hist, v_expected, rtol=hist_tol_all)
+ hist_tol_center = hist_tol_all/10
remove = 3
np.testing.assert_allclose(
- v_hist[remove:-remove], v_expected[remove:-remove], rtol=0.01
+ v_hist[remove:-remove], v_expected[remove:-remove], rtol=hist_tol_center
)
# pressure tensor against expressions from
@@ -291,9 +293,10 @@ def test_resting_fluid(zero_centered: bool, domain_size: int, target=Target.CPU)
# thermalized, the expectation value of <u,u> = kT due to the
# equi-partition theorem.
p_av_expected = np.diag([rho_0 * c_s**2 + kT] * 3)
+ pressure_atol = c_s**2 / 200 *(16/domain_size)**(3/2)
np.testing.assert_allclose(
- np.mean(res_pressure, axis=0), p_av_expected, atol=c_s**2 / 2000
- )
+ np.mean(res_pressure, axis=0), p_av_expected, atol=pressure_atol)
+
@pytest.mark.parametrize(
@@ -314,8 +317,8 @@ def test_point_force(zero_centered: bool, domain_size: int, target=Target.CPU):
kT=kT,
omega_shear=0.8,
omega_bulk=0.5,
- omega_even=0.04,
- omega_odd=0.3,
+ omega_even=0.8,
+ omega_odd=0.8,
zero_centered=zero_centered
)
@@ -326,7 +329,7 @@ def test_point_force(zero_centered: bool, domain_size: int, target=Target.CPU):
introduced_momentum = np.zeros(3)
for i in range(100):
- point_force = 1e-5 * (np.random.random(3) - 0.5)
+ point_force = 1e-2/domain_size**(3/2) * (np.random.random(3) - 0.5)
introduced_momentum += point_force
# Note that ghost layers are included in the indexing
@@ -336,9 +339,10 @@ def test_point_force(zero_centered: bool, domain_size: int, target=Target.CPU):
t = time_loop(t, 1)
res_u = dh.gather_array("u").reshape((-1, 3))
res_rho = dh.gather_array("rho").reshape((-1,))
+
# mass conservation
- np.testing.assert_allclose(np.mean(res_rho), rho_0, atol=3e-12)
+ np.testing.assert_allclose(np.mean(res_rho), rho_0, atol=3E-12)
# momentum conservation
momentum = np.dot(res_rho, res_u)
--
GitLab