Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Commits on Source (6)
Hide whitespace changes
Inline Side-by-side
.gitlab-ci.yml
View file @ 32380aab
...@@ -29,7 +29,7 @@ stages: ...@@ -29,7 +29,7 @@ stages:
tests-and-coverage: tests-and-coverage:
stage: pretest stage: pretest
extends: .every-commit extends: .every-commit
image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full:cupy12.3
script: script:
# - pip install sympy --upgrade # - pip install sympy --upgrade
- export NUM_CORES=$(nproc --all) - export NUM_CORES=$(nproc --all)
......
COPYING.txt
View file @ 32380aab
------------------------ Important ---------------------------------
lbmpy is under the following GNU AGPLv3 license.
This license holds for the sources of lbmpy itself as well
as for all kernels generated with lbmpy i.e.
the output of lbmpy is also protected by the GNU AGPLv3 license.
----------------------------------------------------------------------
GNU AFFERO GENERAL PUBLIC LICENSE GNU AFFERO GENERAL PUBLIC LICENSE
Version 3, 19 November 2007 Version 3, 19 November 2007
......
src/lbmpy/boundaries/wall_function_models.py
View file @ 32380aab
...@@ -157,7 +157,7 @@ class MuskerLaw(ImplicitWallFunctionModel): ...@@ -157,7 +157,7 @@ class MuskerLaw(ImplicitWallFunctionModel):
def law(u_p, y_p): def law(u_p, y_p):
arctan = sp.Float(5.424) * sp.atan(sp.Float(0.119760479041916168) * y_p - sp.Float(0.488023952095808383)) arctan = sp.Float(5.424) * sp.atan(sp.Float(0.119760479041916168) * y_p - sp.Float(0.488023952095808383))
logarithm = (sp.Float(0.434) * sp.log((y_p + sp.Float(10.6)) ** sp.Float(9.6) logarithm = (sp.Float(0.434) * sp.log((y_p + sp.Float(10.6)) ** sp.Float(9.6)
/ (y_p ** 2 - sp.Float(8.15) * y_p + sp.Float(86)) ** 2, 10)) / (y_p ** 2 - sp.Float(8.15) * y_p + sp.Float(86)) ** 2))
return (arctan + logarithm - sp.Float(3.50727901936264842)) - u_p return (arctan + logarithm - sp.Float(3.50727901936264842)) - u_p
u_plus = velocity_symbol / self.u_tau[0] u_plus = velocity_symbol / self.u_tau[0]
......
src/lbmpy/creationfunctions.py
View file @ 32380aab
...@@ -376,8 +376,8 @@ class LBMConfig: ...@@ -376,8 +376,8 @@ class LBMConfig:
if not self.compressible and self.method in (Method.MONOMIAL_CUMULANT, Method.CUMULANT): if not self.compressible and self.method in (Method.MONOMIAL_CUMULANT, Method.CUMULANT):
raise ValueError("Incompressible cumulant-based methods are not supported (yet).") raise ValueError("Incompressible cumulant-based methods are not supported (yet).")
if self.zero_centered and (self.entropic or self.fluctuating): if self.zero_centered and self.entropic:
raise ValueError("Entropic and fluctuating methods can only be created with `zero_centered=False`.") raise ValueError("Entropic methods can only be created with `zero_centered=False`.")
# Check or infer delta-equilibrium # Check or infer delta-equilibrium
if self.delta_equilibrium is not None: if self.delta_equilibrium is not None:
......
src/lbmpy/fluctuatinglb.py
View file @ 32380aab
...@@ -19,9 +19,7 @@ def add_fluctuations_to_collision_rule(collision_rule, temperature=None, amplitu ...@@ -19,9 +19,7 @@ def add_fluctuations_to_collision_rule(collision_rule, temperature=None, amplitu
"""""" """"""
if not (temperature and not amplitudes) or (temperature and amplitudes): if not (temperature and not amplitudes) or (temperature and amplitudes):
raise ValueError("Fluctuating LBM: Pass either 'temperature' or 'amplitudes'.") raise ValueError("Fluctuating LBM: Pass either 'temperature' or 'amplitudes'.")
if collision_rule.method.conserved_quantity_computation.zero_centered_pdfs:
raise ValueError("The fluctuating LBM is not implemented for zero-centered PDF storage.")
method = collision_rule.method method = collision_rule.method
if not amplitudes: if not amplitudes:
amplitudes = fluctuation_amplitude_from_temperature(method, temperature, c_s_sq) amplitudes = fluctuation_amplitude_from_temperature(method, temperature, c_s_sq)
...@@ -44,9 +42,7 @@ def fluctuation_amplitude_from_temperature(method, temperature, c_s_sq=sp.Symbol ...@@ -44,9 +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""" """Produces amplitude equations according to (2.60) and (3.54) in Schiller08"""
normalization_factors = sp.matrix_multiply_elementwise(method.moment_matrix, method.moment_matrix) * \ normalization_factors = sp.matrix_multiply_elementwise(method.moment_matrix, method.moment_matrix) * \
sp.Matrix(method.weights) sp.Matrix(method.weights)
density = method.zeroth_order_equilibrium_moment_symbol density = method._cqc.density_symbol
if method.conserved_quantity_computation.zero_centered_pdfs:
density += 1
mu = temperature * density / c_s_sq mu = temperature * density / c_s_sq
return [sp.sqrt(mu * norm * (1 - (1 - rr) ** 2)) return [sp.sqrt(mu * norm * (1 - (1 - rr) ** 2))
for norm, rr in zip(normalization_factors, method.relaxation_rates)] for norm, rr in zip(normalization_factors, method.relaxation_rates)]
......
tests/test_fluctuating_lb.py
View file @ 32380aab
...@@ -5,13 +5,16 @@ import pytest ...@@ -5,13 +5,16 @@ import pytest
import pystencils as ps import pystencils as ps
from pystencils import get_code_str from pystencils import get_code_str
from pystencils.backends.simd_instruction_sets import get_supported_instruction_sets, get_vector_instruction_set from pystencils.backends.simd_instruction_sets import (
get_supported_instruction_sets,
get_vector_instruction_set,
)
from pystencils.cpu.cpujit import get_compiler_config from pystencils.cpu.cpujit import get_compiler_config
from pystencils.enums import Target from pystencils.enums import Target
from pystencils.rng import PhiloxTwoDoubles from pystencils.rng import PhiloxTwoDoubles
from lbmpy.creationfunctions import * from lbmpy.creationfunctions import *
from lbmpy.forcemodels import Guo from lbmpy.forcemodels import Guo
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
import numpy as np import numpy as np
from lbmpy.enums import Stencil from lbmpy.enums import Stencil
...@@ -20,13 +23,18 @@ from lbmpy.stencils import LBStencil ...@@ -20,13 +23,18 @@ from lbmpy.stencils import LBStencil
def _skip_instruction_sets_windows(instruction_sets): def _skip_instruction_sets_windows(instruction_sets):
if get_compiler_config()['os'] == 'windows': if get_compiler_config()["os"] == "windows":
# skip instruction sets supported by the CPU but not by the compiler # skip instruction sets supported by the CPU but not by the compiler
if 'avx' in instruction_sets and ('/arch:avx2' not in get_compiler_config()['flags'].lower() if "avx" in instruction_sets and (
and '/arch:avx512' not in get_compiler_config()['flags'].lower()): "/arch:avx2" not in get_compiler_config()["flags"].lower()
instruction_sets.remove('avx') and "/arch:avx512" not in get_compiler_config()["flags"].lower()
if 'avx512' in instruction_sets and '/arch:avx512' not in get_compiler_config()['flags'].lower(): ):
instruction_sets.remove('avx512') instruction_sets.remove("avx")
if (
"avx512" in instruction_sets
and "/arch:avx512" not in get_compiler_config()["flags"].lower()
):
instruction_sets.remove("avx512")
return instruction_sets return instruction_sets
...@@ -35,11 +43,13 @@ def single_component_maxwell(x1, x2, kT, mass): ...@@ -35,11 +43,13 @@ def single_component_maxwell(x1, x2, kT, mass):
x = np.linspace(x1, x2, 1000) x = np.linspace(x1, x2, 1000)
try: try:
trapezoid = np.trapezoid # since numpy 2.0 trapezoid = np.trapezoid # since numpy 2.0
except AttributeError: except AttributeError:
trapezoid = np.trapz trapezoid = np.trapz
return trapezoid(np.exp(-mass * x ** 2 / (2. * kT)), x) / np.sqrt(2. * np.pi * kT / mass) return trapezoid(np.exp(-mass * x**2 / (2.0 * kT)), x) / np.sqrt(
2.0 * np.pi * kT / mass
)
def rr_getter(moment_group): def rr_getter(moment_group):
...@@ -71,53 +81,86 @@ def second_order_moment_tensor_assignments(function_values, stencil, output_fiel ...@@ -71,53 +81,86 @@ def second_order_moment_tensor_assignments(function_values, stencil, output_fiel
"""Assignments for calculating the pressure tensor""" """Assignments for calculating the pressure tensor"""
assert len(function_values) == len(stencil) assert len(function_values) == len(stencil)
dim = len(stencil[0]) dim = len(stencil[0])
return [ps.Assignment(output_field(i, j), return [
sum(c[i] * c[j] * f for f, c in zip(function_values, stencil))) ps.Assignment(
for i in range(dim) for j in range(dim)] output_field(i, j),
sum(c[i] * c[j] * f for f, c in zip(function_values, stencil)),
)
for i in range(dim)
for j in range(dim)
]
def add_pressure_output_to_collision_rule(collision_rule, pressure_field): def add_pressure_output_to_collision_rule(collision_rule, pressure_field):
pressure_ouput = second_order_moment_tensor_assignments(collision_rule.method.pre_collision_pdf_symbols, pressure_ouput = second_order_moment_tensor_assignments(
collision_rule.method.stencil, pressure_field) collision_rule.method.pre_collision_pdf_symbols,
collision_rule.method.stencil,
pressure_field,
)
collision_rule.main_assignments = collision_rule.main_assignments + pressure_ouput collision_rule.main_assignments = collision_rule.main_assignments + pressure_ouput
def get_fluctuating_lb(size=None, kT=None, def get_fluctuating_lb(
omega_shear=None, omega_bulk=None, omega_odd=None, omega_even=None, size=None,
rho_0=None, target=None): kT=None,
omega_shear=None,
omega_bulk=None,
omega_odd=None,
omega_even=None,
rho_0=None,
target=None,
zero_centered: bool = False,
):
# Parameters # Parameters
stencil = LBStencil(Stencil.D3Q19) stencil = LBStencil(Stencil.D3Q19)
# Setup data handling # Setup data handling
dh = ps.create_data_handling((size,) * stencil.D, periodicity=True, default_target=target) dh = ps.create_data_handling(
src = dh.add_array('src', values_per_cell=stencil.Q, layout='f') (size,) * stencil.D, periodicity=True, default_target=target
dst = dh.add_array_like('dst', 'src') )
rho = dh.add_array('rho', layout='f', latex_name='\\rho', values_per_cell=1) src = dh.add_array("src", values_per_cell=stencil.Q, layout="f")
u = dh.add_array('u', values_per_cell=dh.dim, layout='f') dst = dh.add_array_like("dst", "src")
pressure_field = dh.add_array('pressure', values_per_cell=( rho = dh.add_array("rho", layout="f", latex_name="\\rho", values_per_cell=1)
3, 3), layout='f', gpu=target == Target.GPU) u = dh.add_array("u", values_per_cell=dh.dim, layout="f")
pressure_field = dh.add_array(
"pressure", values_per_cell=(3, 3), layout="f", gpu=target == Target.GPU
)
force_field = dh.add_array( force_field = dh.add_array(
'force', values_per_cell=stencil.D, layout='f', gpu=target == Target.GPU) "force", values_per_cell=stencil.D, layout="f", gpu=target == Target.GPU
)
# Method setup # Method setup
lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, compressible=True, lbm_config = LBMConfig(
weighted=True, zero_centered=False, relaxation_rates=rr_getter, stencil=stencil,
force_model=Guo(force=force_field.center_vector), method=Method.MRT,
fluctuating={'temperature': kT}, compressible=True,
kernel_type='collide_only') weighted=True,
zero_centered=zero_centered,
relaxation_rates=rr_getter,
force_model=Guo(force=force_field.center_vector),
fluctuating={"temperature": kT},
kernel_type="collide_only",
)
lb_method = create_lb_method(lbm_config=lbm_config) lb_method = create_lb_method(lbm_config=lbm_config)
lbm_config.lb_method = lb_method lbm_config.lb_method = lb_method
lbm_opt = LBMOptimisation(symbolic_field=src, cse_global=True)
collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
add_pressure_output_to_collision_rule(collision_rule, pressure_field)
collision = create_lb_update_rule(collision_rule=collision_rule, lbm_opt = LBMOptimisation(symbolic_field=src, cse_global=True)
lbm_config=lbm_config, lbm_optimisation=lbm_opt) collision_rule = create_lb_collision_rule(
stream = create_stream_pull_with_output_kernel(collision.method, src, dst, lbm_config=lbm_config, lbm_optimisation=lbm_opt
{'density': rho, 'velocity': u}) )
# add_pressure_output_to_collision_rule(collision_rule, pressure_field)
collision = create_lb_update_rule(
collision_rule=collision_rule, lbm_config=lbm_config, lbm_optimisation=lbm_opt
)
stream = create_stream_pull_with_output_kernel(
collision.method,
src,
dst,
{"density": rho, "velocity": u, "moment2": pressure_field},
)
config = ps.CreateKernelConfig(cpu_openmp=False, target=dh.default_target) config = ps.CreateKernelConfig(cpu_openmp=False, target=dh.default_target)
...@@ -128,15 +171,18 @@ def get_fluctuating_lb(size=None, kT=None, ...@@ -128,15 +171,18 @@ def get_fluctuating_lb(size=None, kT=None,
sync_pdfs = dh.synchronization_function([src.name]) sync_pdfs = dh.synchronization_function([src.name])
# Initialization # Initialization
init = macroscopic_values_setter(collision.method, velocity=(0,) * dh.dim, init = macroscopic_values_setter(
pdfs=src.center_vector, density=rho.center) collision.method,
velocity=(0,) * dh.dim,
pdfs=src.center_vector,
density=rho_0
)
init_kernel = ps.create_kernel(init, ghost_layers=0).compile() init_kernel = ps.create_kernel(init, ghost_layers=0).compile()
dh.fill(rho.name, rho_0) dh.fill(rho.name, rho_0)
dh.fill(u.name, np.nan, ghost_layers=True, inner_ghost_layers=True) dh.fill(u.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
dh.fill(u.name, 0) dh.fill(u.name, 0)
dh.fill(force_field.name, np.nan, dh.fill(force_field.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
ghost_layers=True, inner_ghost_layers=True)
dh.fill(force_field.name, 0) dh.fill(force_field.name, 0)
dh.run_kernel(init_kernel) dh.run_kernel(init_kernel)
...@@ -144,8 +190,15 @@ def get_fluctuating_lb(size=None, kT=None, ...@@ -144,8 +190,15 @@ def get_fluctuating_lb(size=None, kT=None,
def time_loop(start, steps): def time_loop(start, steps):
dh.all_to_gpu() dh.all_to_gpu()
for i in range(start, start + steps): for i in range(start, start + steps):
dh.run_kernel(collision_kernel, omega_shear=omega_shear, omega_bulk=omega_bulk, dh.run_kernel(
omega_odd=omega_odd, omega_even=omega_even, seed=42, time_step=i) collision_kernel,
omega_shear=omega_shear,
omega_bulk=omega_bulk,
omega_odd=omega_odd,
omega_even=omega_even,
seed=42,
time_step=i,
)
sync_pdfs() sync_pdfs()
dh.run_kernel(stream_kernel) dh.run_kernel(stream_kernel)
...@@ -156,13 +209,27 @@ def get_fluctuating_lb(size=None, kT=None, ...@@ -156,13 +209,27 @@ def get_fluctuating_lb(size=None, kT=None,
return dh, time_loop return dh, time_loop
def test_resting_fluid(target=Target.CPU): @pytest.mark.parametrize(
rho_0 = 0.86 "zero_centered", [False, True], ids=["regular-storage", "zero-centered"]
kT = 4E-4 )
L = [60] * 3 @pytest.mark.parametrize(
dh, time_loop = get_fluctuating_lb(size=L[0], target=target, "domain_size", [8, 60]
rho_0=rho_0, kT=kT, )
omega_shear=0.8, omega_bulk=0.5, omega_even=.04, omega_odd=0.3) def test_resting_fluid(zero_centered: bool, domain_size: int, target=Target.CPU):
rho_0 = 0.86
kT = 4e-4
L = [domain_size] * 3
dh, time_loop = get_fluctuating_lb(
size=L[0],
target=target,
rho_0=rho_0,
kT=kT,
omega_shear=0.8,
omega_bulk=0.5,
omega_even=0.04,
omega_odd=0.3,
zero_centered=zero_centered,
)
# Test # Test
t = 0 t = 0
...@@ -176,38 +243,43 @@ def test_resting_fluid(target=Target.CPU): ...@@ -176,38 +243,43 @@ def test_resting_fluid(target=Target.CPU):
res_u = dh.gather_array("u").reshape((-1, 3)) res_u = dh.gather_array("u").reshape((-1, 3))
res_rho = dh.gather_array("rho").reshape((-1,)) res_rho = dh.gather_array("rho").reshape((-1,))
# mass conservation # mass conservationo
# density per cell fluctuates, but toal mass is conserved
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 conservation
momentum = np.dot(res_rho, res_u) momentum = np.dot(res_rho, res_u)
np.testing.assert_allclose(momentum, [0, 0, 0], atol=1E-10) 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) kinetic_energy = 1 / 2 * np.dot(res_rho, res_u * res_u) / np.prod(L)
np.testing.assert_allclose( kT_tol = 0.075 *(16/domain_size)**(3/2)
kinetic_energy, [kT / 2] * 3, atol=kT * 0.01) np.testing.assert_allclose(kinetic_energy, [kT / 2] * 3, rtol=kT_tol)
# velocity distribution # velocity distribution
v_hist, v_bins = np.histogram( v_hist, v_bins = np.histogram(
res_u, bins=11, range=(-.075, .075), density=True) res_u, bins=11, range=(-0.075, 0.075), density=True
)
# Calculate expected values from single # Calculate expected values from single
v_expected = [] v_expected = []
for j in range(len(v_hist)): for j in range(len(v_hist)):
# Maxwell distribution # Maxwell distribution
res = 1. / (v_bins[j + 1] - v_bins[j]) * \ res = (
single_component_maxwell( 1.0
v_bins[j], v_bins[j + 1], kT, rho_0) / (v_bins[j + 1] - v_bins[j])
* single_component_maxwell(v_bins[j], v_bins[j + 1], kT, rho_0)
)
v_expected.append(res) v_expected.append(res)
v_expected = np.array(v_expected) v_expected = np.array(v_expected)
# 10% accuracy on the entire histogram hist_tol_all = 0.75 *(16/domain_size)**(3/2)
np.testing.assert_allclose(v_hist, v_expected, rtol=0.1) np.testing.assert_allclose(v_hist, v_expected, rtol=hist_tol_all)
# 1% accuracy on the middle part hist_tol_center = hist_tol_all/10
remove = 3 remove = 3
np.testing.assert_allclose( 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 # pressure tensor against expressions from
# Duenweg, Schiller, Ladd, https://arxiv.org/abs/0707.1581 # Duenweg, Schiller, Ladd, https://arxiv.org/abs/0707.1581
...@@ -220,19 +292,35 @@ def test_resting_fluid(target=Target.CPU): ...@@ -220,19 +292,35 @@ def test_resting_fluid(target=Target.CPU):
# Diagonal elements are rho c_s^22 +<u,u>. When the fluid is # Diagonal elements are rho c_s^22 +<u,u>. When the fluid is
# thermalized, the expectation value of <u,u> = kT due to the # thermalized, the expectation value of <u,u> = kT due to the
# equi-partition theorem. # equi-partition theorem.
p_av_expected = np.diag([rho_0 * c_s ** 2 + kT] * 3) 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.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)
def test_point_force(target=Target.CPU):
@pytest.mark.parametrize(
"zero_centered", [False, True], ids=["regular-storage", "zero-centered"]
)
@pytest.mark.parametrize(
"domain_size", [8, 60]
)
def test_point_force(zero_centered: bool, domain_size: int, target=Target.CPU):
"""Test momentum balance for thermalized fluid with applied poitn forces""" """Test momentum balance for thermalized fluid with applied poitn forces"""
rho_0 = 0.86 rho_0 = 0.86
kT = 4E-4 kT = 4e-4
L = [8] * 3 L = [domain_size] * 3
dh, time_loop = get_fluctuating_lb(size=L[0], target=target, dh, time_loop = get_fluctuating_lb(
rho_0=rho_0, kT=kT, size=L[0],
omega_shear=0.8, omega_bulk=0.5, omega_even=.04, omega_odd=0.3) target=target,
rho_0=rho_0,
kT=kT,
omega_shear=0.8,
omega_bulk=0.5,
omega_even=0.8,
omega_odd=0.8,
zero_centered=zero_centered
)
# Test # Test
t = 0 t = 0
...@@ -241,17 +329,17 @@ def test_point_force(target=Target.CPU): ...@@ -241,17 +329,17 @@ def test_point_force(target=Target.CPU):
introduced_momentum = np.zeros(3) introduced_momentum = np.zeros(3)
for i in range(100): for i in range(100):
point_force = 1E-5 * (np.random.random(3) - .5) point_force = 1e-2/domain_size**(3/2) * (np.random.random(3) - 0.5)
introduced_momentum += point_force introduced_momentum += point_force
# Note that ghost layers are included in the indexing # Note that ghost layers are included in the indexing
force_pos = np.random.randint(1, L[0] - 2, size=3) force_pos = np.random.randint(1, L[0] - 2, size=3)
dh.cpu_arrays["force"][force_pos[0], dh.cpu_arrays["force"][force_pos[0], force_pos[1], force_pos[2]] = point_force
force_pos[1], force_pos[2]] = point_force
t = time_loop(t, 1) t = time_loop(t, 1)
res_u = dh.gather_array("u").reshape((-1, 3)) res_u = dh.gather_array("u").reshape((-1, 3))
res_rho = dh.gather_array("rho").reshape((-1,)) res_rho = dh.gather_array("rho").reshape((-1,))
# mass conservation # 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)
...@@ -259,52 +347,72 @@ def test_point_force(target=Target.CPU): ...@@ -259,52 +347,72 @@ def test_point_force(target=Target.CPU):
# momentum conservation # momentum conservation
momentum = np.dot(res_rho, res_u) momentum = np.dot(res_rho, res_u)
np.testing.assert_allclose( np.testing.assert_allclose(
momentum, introduced_momentum + 0.5 * point_force, atol=1E-10) momentum, introduced_momentum + 0.5 * point_force, atol=1e-10
dh.cpu_arrays["force"][force_pos[0], )
force_pos[1], force_pos[2]] = np.zeros(3) dh.cpu_arrays["force"][force_pos[0], force_pos[1], force_pos[2]] = np.zeros(3)
@pytest.mark.skipif(not get_supported_instruction_sets(), reason="No vector instruction sets supported") @pytest.mark.skipif(
@pytest.mark.parametrize('data_type', ("float32", "float64")) not get_supported_instruction_sets(), reason="No vector instruction sets supported"
@pytest.mark.parametrize('assume_aligned', (True, False)) )
@pytest.mark.parametrize('assume_inner_stride_one', (True, False)) @pytest.mark.parametrize("data_type", ("float32", "float64"))
@pytest.mark.parametrize('assume_sufficient_line_padding', (True, False)) @pytest.mark.parametrize("assume_aligned", (True, False))
def test_vectorization(data_type, assume_aligned, assume_inner_stride_one, assume_sufficient_line_padding): @pytest.mark.parametrize("assume_inner_stride_one", (True, False))
@pytest.mark.parametrize("assume_sufficient_line_padding", (True, False))
def test_vectorization(
data_type, assume_aligned, assume_inner_stride_one, assume_sufficient_line_padding
):
stencil = LBStencil(Stencil.D3Q19) stencil = LBStencil(Stencil.D3Q19)
pdfs, pdfs_tmp = ps.fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[3D]", layout='fzyx') pdfs, pdfs_tmp = ps.fields(
f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[3D]", layout="fzyx"
)
method = create_mrt_orthogonal( method = create_mrt_orthogonal(
stencil=stencil, stencil=stencil, compressible=True, weighted=True, relaxation_rates=rr_getter
)
rng_node = (
ps.rng.PhiloxTwoDoubles if data_type == "float64" else ps.rng.PhiloxFourFloats
)
lbm_config = LBMConfig(
lb_method=method,
fluctuating={
"temperature": sp.Symbol("kT"),
"rng_node": rng_node,
"block_offsets": tuple([0] * stencil.D),
},
compressible=True, compressible=True,
weighted=True, zero_centered=False,
relaxation_rates=rr_getter) stencil=method.stencil,
kernel_type="collide_only",
rng_node = ps.rng.PhiloxTwoDoubles if data_type == "float64" else ps.rng.PhiloxFourFloats )
lbm_config = LBMConfig(lb_method=method, fluctuating={'temperature': sp.Symbol("kT"), lbm_opt = LBMOptimisation(
'rng_node': rng_node, cse_global=True, symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp
'block_offsets': tuple([0] * stencil.D)}, )
compressible=True, zero_centered=False,
stencil=method.stencil, kernel_type='collide_only')
lbm_opt = LBMOptimisation(cse_global=True, symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp)
collision = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt) collision = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
instruction_sets = _skip_instruction_sets_windows(get_supported_instruction_sets()) instruction_sets = _skip_instruction_sets_windows(get_supported_instruction_sets())
instruction_set = instruction_sets[-1] instruction_set = instruction_sets[-1]
config = ps.CreateKernelConfig(target=Target.CPU, config = ps.CreateKernelConfig(
data_type=data_type, default_number_float=data_type, target=Target.CPU,
cpu_vectorize_info={'instruction_set': instruction_set, data_type=data_type,
'assume_aligned': assume_aligned, default_number_float=data_type,
'assume_inner_stride_one': assume_inner_stride_one, cpu_vectorize_info={
'assume_sufficient_line_padding': assume_sufficient_line_padding, "instruction_set": instruction_set,
} "assume_aligned": assume_aligned,
) "assume_inner_stride_one": assume_inner_stride_one,
"assume_sufficient_line_padding": assume_sufficient_line_padding,
if not assume_inner_stride_one and 'storeS' not in get_vector_instruction_set(data_type, instruction_set): },
)
if not assume_inner_stride_one and "storeS" not in get_vector_instruction_set(
data_type, instruction_set
):
with pytest.warns(UserWarning) as pytest_warnings: with pytest.warns(UserWarning) as pytest_warnings:
ast = ps.create_kernel(collision, config=config) ast = ps.create_kernel(collision, config=config)
assert 'Could not vectorize loop' in pytest_warnings[0].message.args[0] assert "Could not vectorize loop" in pytest_warnings[0].message.args[0]
else: else:
ast = ps.create_kernel(collision, config=config) ast = ps.create_kernel(collision, config=config)
ast.compile() ast.compile()
...@@ -312,31 +420,54 @@ def test_vectorization(data_type, assume_aligned, assume_inner_stride_one, assum ...@@ -312,31 +420,54 @@ def test_vectorization(data_type, assume_aligned, assume_inner_stride_one, assum
print(code) print(code)
@pytest.mark.parametrize('data_type', ("float32", "float64")) @pytest.mark.parametrize("data_type", ("float32", "float64"))
@pytest.mark.parametrize('assume_aligned', (True, False)) @pytest.mark.parametrize("assume_aligned", (True, False))
@pytest.mark.parametrize('assume_inner_stride_one', (True, False)) @pytest.mark.parametrize("assume_inner_stride_one", (True, False))
@pytest.mark.parametrize('assume_sufficient_line_padding', (True, False)) @pytest.mark.parametrize("assume_sufficient_line_padding", (True, False))
def test_fluctuating_lb_issue_188_wlb(data_type, assume_aligned, def test_fluctuating_lb_issue_188_wlb(
assume_inner_stride_one, assume_sufficient_line_padding): data_type, assume_aligned, assume_inner_stride_one, assume_sufficient_line_padding
):
stencil = LBStencil(Stencil.D3Q19) stencil = LBStencil(Stencil.D3Q19)
temperature = sp.symbols("temperature") temperature = sp.symbols("temperature")
pdfs, pdfs_tmp = ps.fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[3D]", layout='fzyx') pdfs, pdfs_tmp = ps.fields(
f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[3D]", layout="fzyx"
rng_node = ps.rng.PhiloxTwoDoubles if data_type == "float64" else ps.rng.PhiloxFourFloats )
fluctuating = {'temperature': temperature,
'block_offsets': 'walberla', rng_node = (
'rng_node': rng_node} ps.rng.PhiloxTwoDoubles if data_type == "float64" else ps.rng.PhiloxFourFloats
)
lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, compressible=True, fluctuating = {
weighted=True, zero_centered=False, relaxation_rate=1.4, "temperature": temperature,
fluctuating=fluctuating) "block_offsets": "walberla",
lbm_opt = LBMOptimisation(symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp, cse_global=True) "rng_node": rng_node,
}
lbm_config = LBMConfig(
stencil=stencil,
method=Method.MRT,
compressible=True,
weighted=True,
zero_centered=False,
relaxation_rate=1.4,
fluctuating=fluctuating,
)
lbm_opt = LBMOptimisation(
symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp, cse_global=True
)
up = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt) up = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
cpu_vectorize_info = {'instruction_set': 'avx', 'assume_inner_stride_one': True, 'assume_aligned': True} cpu_vectorize_info = {
config = ps.CreateKernelConfig(target=ps.Target.CPU, data_type=data_type, default_number_float=data_type, "instruction_set": "avx",
cpu_vectorize_info=cpu_vectorize_info) "assume_inner_stride_one": True,
"assume_aligned": True,
}
config = ps.CreateKernelConfig(
target=ps.Target.CPU,
data_type=data_type,
default_number_float=data_type,
cpu_vectorize_info=cpu_vectorize_info,
)
ast = create_kernel(up, config=config) ast = create_kernel(up, config=config)
code = ps.get_code_str(ast) code = ps.get_code_str(ast)
......