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
Select Git revision
  • Sparse
  • WallLaw
  • improved_comm
  • master
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
  • release/1.3.7
42 results

Target

Select target project
  • ravi.k.ayyala/lbmpy
  • brendan-waters/lbmpy
  • anirudh.jonnalagadda/lbmpy
  • jbadwaik/lbmpy
  • alexander.reinauer/lbmpy
  • itischler/lbmpy
  • he66coqe/lbmpy
  • ev81oxyl/lbmpy
  • Bindgen/lbmpy
  • da15siwa/lbmpy
  • holzer/lbmpy
  • RudolfWeeber/lbmpy
  • pycodegen/lbmpy
13 results
Select Git revision
  • GetterSetterAPI
  • HRR
  • HydroPressure
  • InplaceConfig
  • Outflow
  • PhaseField
  • Sparse
  • UBBVelocity
  • UpdateAPISparse
  • WallLaw
  • WetNodeBoundaries
  • csebug
  • feature/sparse
  • feature/try
  • improved_comm
  • install_requires
  • master
  • phaseField
  • relaxationrates
  • test_martin
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
57 results
Show changes
Showing
with 357 additions and 116 deletions
...@@ -135,8 +135,7 @@ def get_triple_points(phase_arr, phase_indices, contour_line_eps=0.01, threshold ...@@ -135,8 +135,7 @@ def get_triple_points(phase_arr, phase_indices, contour_line_eps=0.01, threshold
def analytic_neumann_angles(kappas): def analytic_neumann_angles(kappas):
"""Computes analytic Neumann angles using surface tension parameters. """Computes analytic Neumann angles using surface tension parameters.
>>> analytic_neumann_angles([0.1, 0.1, 0.1]) >>> assert analytic_neumann_angles([0.1, 0.1, 0.1]) == [120.00000000000001, 120.00000000000001, 120.00000000000001]
[120.00000000000001, 120.00000000000001, 120.00000000000001]
>>> r = analytic_neumann_angles([0.1, 0.2, 0.3]) >>> r = analytic_neumann_angles([0.1, 0.2, 0.3])
>>> assert np.allclose(sum(r), 360) >>> assert np.allclose(sum(r), 360)
......
import math import math
import sympy as sp import sympy as sp
from .._compat import IS_PYSTENCILS_2
if IS_PYSTENCILS_2:
raise ImportError("`lbmpy.phasefield_allen_cahn.contact_angle` is only available when running with pystencils 1.x")
from pystencils.astnodes import Block, Conditional, SympyAssignment from pystencils.astnodes import Block, Conditional, SympyAssignment
from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
from pystencils.boundaries.boundaryconditions import Boundary from pystencils.boundaries.boundaryconditions import Boundary
from pystencils.typing import TypedSymbol from pystencils import TypedSymbol
from pystencils.typing import CastFunc from pystencils.typing import CastFunc
......
from pystencils import Assignment, AssignmentCollection from pystencils import Assignment, AssignmentCollection
from pystencils.sympyextensions import scalar_product from pystencils.sympyextensions import scalar_product
from pystencils.simp.subexpression_insertion import insert_constants from pystencils.simp import insert_constants
from lbmpy.phasefield_allen_cahn.derivatives import isotropic_gradient_symbolic, laplacian_symbolic from lbmpy.phasefield_allen_cahn.derivatives import isotropic_gradient_symbolic, laplacian_symbolic
......
...@@ -366,8 +366,8 @@ def calculate_parameters_taylor_bubble(reference_length=128, ...@@ -366,8 +366,8 @@ def calculate_parameters_taylor_bubble(reference_length=128,
water_mu_lattice_units = math.sqrt(density_diff * density_heavy * grav_df_cube) / inverse_viscosity_number water_mu_lattice_units = math.sqrt(density_diff * density_heavy * grav_df_cube) / inverse_viscosity_number
air_mu_lattice_units = water_mu_lattice_units / (water_mu / air_mu) air_mu_lattice_units = water_mu_lattice_units / (water_mu / air_mu)
dynamic_viscosity_heavy = water_mu_lattice_units / density_heavy dynamic_viscosity_heavy = water_mu_lattice_units
dynamic_viscosity_light = air_mu_lattice_units / density_light dynamic_viscosity_light = air_mu_lattice_units
surface_tension_lattice_units = density_diff * gravity_lattice_units * diameter_fluid ** 2 / bond_number surface_tension_lattice_units = density_diff * gravity_lattice_units * diameter_fluid ** 2 / bond_number
......
...@@ -27,9 +27,11 @@ that defines the viscosity of the fluid (valid values being between 0 and 2). ...@@ -27,9 +27,11 @@ that defines the viscosity of the fluid (valid values being between 0 and 2).
""" """
import numpy as np import numpy as np
from lbmpy._compat import IS_PYSTENCILS_2
from lbmpy.boundaries import UBB, FixedDensity, NoSlip from lbmpy.boundaries import UBB, FixedDensity, NoSlip
from lbmpy.geometry import add_pipe_inflow_boundary, add_pipe_walls from lbmpy.geometry import add_pipe_inflow_boundary, add_pipe_walls
from lbmpy.lbstep import LatticeBoltzmannStep from lbmpy.lbstep import LatticeBoltzmannStep
from pystencils import Target
from pystencils.datahandling import create_data_handling from pystencils.datahandling import create_data_handling
from pystencils.slicing import slice_from_direction from pystencils.slicing import slice_from_direction
...@@ -86,7 +88,7 @@ def create_lid_driven_cavity(domain_size=None, lid_velocity=0.005, lbm_kernel=No ...@@ -86,7 +88,7 @@ def create_lid_driven_cavity(domain_size=None, lid_velocity=0.005, lbm_kernel=No
assert domain_size is not None or data_handling is not None assert domain_size is not None or data_handling is not None
if data_handling is None: if data_handling is None:
optimization = kwargs.get('optimization', None) optimization = kwargs.get('optimization', None)
target = optimization.get('target', None) if optimization else None target = optimization.get('target', None) if optimization else Target.CPU
data_handling = create_data_handling(domain_size, data_handling = create_data_handling(domain_size,
periodicity=False, periodicity=False,
default_ghost_layers=1, default_ghost_layers=1,
...@@ -94,7 +96,13 @@ def create_lid_driven_cavity(domain_size=None, lid_velocity=0.005, lbm_kernel=No ...@@ -94,7 +96,13 @@ def create_lid_driven_cavity(domain_size=None, lid_velocity=0.005, lbm_kernel=No
default_target=target) default_target=target)
step = LatticeBoltzmannStep(data_handling=data_handling, lbm_kernel=lbm_kernel, name="ldc", **kwargs) step = LatticeBoltzmannStep(data_handling=data_handling, lbm_kernel=lbm_kernel, name="ldc", **kwargs)
my_ubb = UBB(velocity=[lid_velocity, 0, 0][:step.method.dim]) if IS_PYSTENCILS_2:
my_ubb = UBB(
velocity=[lid_velocity, 0, 0][:step.method.dim],
data_type=step.config.get_option("default_dtype")
)
else: # legacy version does not set data type of boundary correctly
my_ubb = UBB(velocity=[lid_velocity, 0, 0][:step.method.dim])
step.boundary_handling.set_boundary(my_ubb, slice_from_direction('N', step.dim)) step.boundary_handling.set_boundary(my_ubb, slice_from_direction('N', step.dim))
for direction in ('W', 'E', 'S') if step.dim == 2 else ('W', 'E', 'S', 'T', 'B'): for direction in ('W', 'E', 'S') if step.dim == 2 else ('W', 'E', 'S', 'T', 'B'):
step.boundary_handling.set_boundary(NoSlip(), slice_from_direction(direction, step.dim)) step.boundary_handling.set_boundary(NoSlip(), slice_from_direction(direction, step.dim))
......
...@@ -3,7 +3,13 @@ from typing import Tuple ...@@ -3,7 +3,13 @@ from typing import Tuple
import numpy as np import numpy as np
import sympy as sp import sympy as sp
from lbmpy.custom_code_nodes import LbmWeightInfo from .._compat import IS_PYSTENCILS_2
if IS_PYSTENCILS_2:
from lbmpy.lookup_tables import LbmWeightInfo
else:
from lbmpy.custom_code_nodes import LbmWeightInfo
from pystencils import Assignment, Field, TypedSymbol from pystencils import Assignment, Field, TypedSymbol
from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
from pystencils.boundaries.createindexlist import ( from pystencils.boundaries.createindexlist import (
...@@ -123,7 +129,10 @@ class SparseLbMapper: ...@@ -123,7 +129,10 @@ class SparseLbMapper:
assert direction_idx == 0 assert direction_idx == 0
continue continue
for own_cell_idx, cell in enumerate(self.fluid_coordinates): for own_cell_idx, cell in enumerate(self.fluid_coordinates):
inv_neighbor_cell = np.array([cell_i - dir_i for cell_i, dir_i in zip(cell, direction)]) inv_neighbor_cell = np.array(
[np.int32(cell_i) - dir_i for cell_i, dir_i in zip(cell, direction)],
dtype=np.uint32
)
if flag_arr[tuple(inv_neighbor_cell)] & fluid_boundary_mask: if flag_arr[tuple(inv_neighbor_cell)] & fluid_boundary_mask:
neighbor_cell_idx = self.cell_idx(tuple(inv_neighbor_cell)) neighbor_cell_idx = self.cell_idx(tuple(inv_neighbor_cell))
result.append(pdf_index(neighbor_cell_idx, direction_idx)) result.append(pdf_index(neighbor_cell_idx, direction_idx))
...@@ -232,7 +241,17 @@ class SparseLbBoundaryMapper: ...@@ -232,7 +241,17 @@ class SparseLbBoundaryMapper:
return result return result
def assignments(self): def assignments(self):
return [BoundaryOffsetInfo(self.method.stencil), if IS_PYSTENCILS_2:
return (
BoundaryOffsetInfo(self.method.stencil).get_array_declarations()
+ LbmWeightInfo(self.method).get_array_declarations()
+ [Assignment(self.DIR_SYMBOL, self.index_field(self.DIR_SYMBOL.name))]
+ self.boundary_eqs
)
else:
return [
BoundaryOffsetInfo(self.method.stencil),
LbmWeightInfo(self.method), LbmWeightInfo(self.method),
Assignment(self.DIR_SYMBOL, self.index_field(self.DIR_SYMBOL.name)), Assignment(self.DIR_SYMBOL, self.index_field(self.DIR_SYMBOL.name)),
*self.boundary_eqs] *self.boundary_eqs
]
...@@ -120,6 +120,26 @@ def create_stream_pull_with_output_kernel(lb_method, src_field, dst_field=None, ...@@ -120,6 +120,26 @@ def create_stream_pull_with_output_kernel(lb_method, src_field, dst_field=None,
simplification_hints=output_eq_collection.simplification_hints) simplification_hints=output_eq_collection.simplification_hints)
def create_copy_kernel(stencil, src_field, dst_field, accessor=StreamPullTwoFieldsAccessor()):
"""Creates a copy kernel, which can be used to transfer information from src to dst field.
Args:
stencil: lattice Boltzmann stencil which is used in the form of a tuple of tuples
src_field: field used for reading pdf values
dst_field: field used for writing pdf values
accessor: instance of PdfFieldAccessor, defining where to read and write values
to create e.g. a fused stream-collide kernel See 'fieldaccess.PdfFieldAccessor'
Returns:
AssignmentCollection of a copy update rule
"""
temporary_symbols = sp.symbols(f'copied:{stencil.Q}')
subexpressions = [Assignment(tmp, acc) for tmp, acc in zip(temporary_symbols, accessor.write(src_field, stencil))]
main_assignments = [Assignment(acc, tmp) for acc, tmp in zip(accessor.write(dst_field, stencil), temporary_symbols)]
return AssignmentCollection(main_assignments, subexpressions=subexpressions)
# ---------------------------------- Pdf array creation for various layouts -------------------------------------------- # ---------------------------------- Pdf array creation for various layouts --------------------------------------------
......
...@@ -11,7 +11,7 @@ from lbmpy.stencils import LBStencil ...@@ -11,7 +11,7 @@ from lbmpy.stencils import LBStencil
import pystencils as ps import pystencils as ps
from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object
from pystencils.typing import TypedSymbol, create_type from pystencils import TypedSymbol
from pystencils.field import Field, FieldType from pystencils.field import Field, FieldType
import pytest import pytest
...@@ -43,7 +43,7 @@ def test_advanced_streaming_noslip_single_cell(stencil, streaming_pattern, prev_ ...@@ -43,7 +43,7 @@ def test_advanced_streaming_noslip_single_cell(stencil, streaming_pattern, prev_
index_struct_dtype = numpy_data_type_for_boundary_object(noslip, stencil.D) index_struct_dtype = numpy_data_type_for_boundary_object(noslip, stencil.D)
index_field = Field('indexVector', FieldType.INDEXED, index_struct_dtype, layout=[0], index_field = Field('indexVector', FieldType.INDEXED, index_struct_dtype, layout=[0],
shape=(TypedSymbol("indexVectorSize", create_type(np.int64)), 1), strides=(1, 1)) shape=(TypedSymbol("indexVectorSize", np.int32), 1), strides=(1, 1))
index_vector = np.array([pos + (d,) for d in range(stencil.Q)], dtype=index_struct_dtype) index_vector = np.array([pos + (d,) for d in range(stencil.Q)], dtype=index_struct_dtype)
ast = create_lattice_boltzmann_boundary_kernel(pdf_field, ast = create_lattice_boltzmann_boundary_kernel(pdf_field,
......
...@@ -5,36 +5,52 @@ import numpy as np ...@@ -5,36 +5,52 @@ import numpy as np
from lbmpy.stencils import LBStencil from lbmpy.stencils import LBStencil
from pystencils.slicing import get_slice_before_ghost_layer, get_ghost_region_slice from pystencils.slicing import get_slice_before_ghost_layer, get_ghost_region_slice
from lbmpy.creationfunctions import create_lb_update_rule, LBMConfig, LBMOptimisation from lbmpy.creationfunctions import create_lb_update_rule, LBMConfig, LBMOptimisation
from lbmpy.advanced_streaming.communication import get_communication_slices, _fix_length_one_slices, \ from lbmpy.advanced_streaming.communication import (
LBMPeriodicityHandling get_communication_slices,
_fix_length_one_slices,
LBMPeriodicityHandling,
periodic_pdf_gpu_copy_kernel,
)
from lbmpy.advanced_streaming.utility import streaming_patterns, Timestep from lbmpy.advanced_streaming.utility import streaming_patterns, Timestep
from lbmpy.enums import Stencil from lbmpy.enums import Stencil
import pytest import pytest
@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]) @pytest.mark.parametrize(
@pytest.mark.parametrize('streaming_pattern', streaming_patterns) "stencil", [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]
@pytest.mark.parametrize('timestep', [Timestep.EVEN, Timestep.ODD]) )
@pytest.mark.parametrize("streaming_pattern", streaming_patterns)
@pytest.mark.parametrize("timestep", [Timestep.EVEN, Timestep.ODD])
def test_slices_not_empty(stencil, streaming_pattern, timestep): def test_slices_not_empty(stencil, streaming_pattern, timestep):
stencil = LBStencil(stencil) stencil = LBStencil(stencil)
arr = np.zeros((4,) * stencil.D + (stencil.Q,)) arr = np.zeros((4,) * stencil.D + (stencil.Q,))
slices = get_communication_slices(stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep, slices = get_communication_slices(
ghost_layers=1) stencil,
streaming_pattern=streaming_pattern,
prev_timestep=timestep,
ghost_layers=1,
)
for _, slices_list in slices.items(): for _, slices_list in slices.items():
for src, dst in slices_list: for src, dst in slices_list:
assert all(s != 0 for s in arr[src].shape) assert all(s != 0 for s in arr[src].shape)
assert all(s != 0 for s in arr[dst].shape) assert all(s != 0 for s in arr[dst].shape)
@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]) @pytest.mark.parametrize(
@pytest.mark.parametrize('streaming_pattern', streaming_patterns) "stencil", [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]
@pytest.mark.parametrize('timestep', [Timestep.EVEN, Timestep.ODD]) )
@pytest.mark.parametrize("streaming_pattern", streaming_patterns)
@pytest.mark.parametrize("timestep", [Timestep.EVEN, Timestep.ODD])
def test_src_dst_same_shape(stencil, streaming_pattern, timestep): def test_src_dst_same_shape(stencil, streaming_pattern, timestep):
stencil = LBStencil(stencil) stencil = LBStencil(stencil)
arr = np.zeros((4,) * stencil.D + (stencil.Q,)) arr = np.zeros((4,) * stencil.D + (stencil.Q,))
slices = get_communication_slices(stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep, slices = get_communication_slices(
ghost_layers=1) stencil,
streaming_pattern=streaming_pattern,
prev_timestep=timestep,
ghost_layers=1,
)
for _, slices_list in slices.items(): for _, slices_list in slices.items():
for src, dst in slices_list: for src, dst in slices_list:
src_shape = arr[src].shape src_shape = arr[src].shape
...@@ -42,12 +58,15 @@ def test_src_dst_same_shape(stencil, streaming_pattern, timestep): ...@@ -42,12 +58,15 @@ def test_src_dst_same_shape(stencil, streaming_pattern, timestep):
assert src_shape == dst_shape assert src_shape == dst_shape
@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]) @pytest.mark.parametrize(
"stencil", [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]
)
def test_pull_communication_slices(stencil): def test_pull_communication_slices(stencil):
stencil = LBStencil(stencil) stencil = LBStencil(stencil)
slices = get_communication_slices( slices = get_communication_slices(
stencil, streaming_pattern='pull', prev_timestep=Timestep.BOTH, ghost_layers=1) stencil, streaming_pattern="pull", prev_timestep=Timestep.BOTH, ghost_layers=1
)
for i, d in enumerate(stencil): for i, d in enumerate(stencil):
if i == 0: if i == 0:
continue continue
...@@ -58,21 +77,115 @@ def test_pull_communication_slices(stencil): ...@@ -58,21 +77,115 @@ def test_pull_communication_slices(stencil):
dst = s[1][:-1] dst = s[1][:-1]
break break
inner_slice = _fix_length_one_slices(get_slice_before_ghost_layer(d, ghost_layers=1)) inner_slice = _fix_length_one_slices(
get_slice_before_ghost_layer(d, ghost_layers=1)
)
inv_dir = (-e for e in d) inv_dir = (-e for e in d)
gl_slice = _fix_length_one_slices(get_ghost_region_slice(inv_dir, ghost_layers=1)) gl_slice = _fix_length_one_slices(
get_ghost_region_slice(inv_dir, ghost_layers=1)
)
assert src == inner_slice assert src == inner_slice
assert dst == gl_slice assert dst == gl_slice
@pytest.mark.parametrize('stencil_name', [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]) @pytest.mark.parametrize("direction", LBStencil(Stencil.D3Q27).stencil_entries)
@pytest.mark.parametrize("pull", [False, True])
def test_gpu_comm_kernels(direction: tuple, pull: bool):
pytest.importorskip("cupy")
stencil = LBStencil(Stencil.D3Q27)
inv_dir = stencil[stencil.inverse_index(direction)]
target = ps.Target.GPU
domain_size = (4,) * stencil.D
dh: ps.datahandling.SerialDataHandling = ps.create_data_handling(
domain_size,
periodicity=(True,) * stencil.D,
parallel=False,
default_target=target,
)
field = dh.add_array("field", values_per_cell=2)
if pull:
dst_slice = get_ghost_region_slice(inv_dir)
src_slice = get_slice_before_ghost_layer(direction)
else:
dst_slice = get_slice_before_ghost_layer(direction)
src_slice = get_ghost_region_slice(inv_dir)
src_slice += (1,)
dst_slice += (1,)
kernel = periodic_pdf_gpu_copy_kernel(field, src_slice, dst_slice)
dh.cpu_arrays[field.name][src_slice] = 42.0
dh.all_to_gpu()
dh.run_kernel(kernel)
dh.all_to_cpu()
np.testing.assert_equal(dh.cpu_arrays[field.name][dst_slice], 42.0)
@pytest.mark.parametrize("stencil", [Stencil.D2Q9, Stencil.D3Q19])
@pytest.mark.parametrize("streaming_pattern", streaming_patterns)
def test_direct_copy_and_kernels_equivalence(stencil: Stencil, streaming_pattern: str):
pytest.importorskip("cupy")
target = ps.Target.GPU
stencil = LBStencil(stencil)
domain_size = (4,) * stencil.D
dh: ps.datahandling.SerialDataHandling = ps.create_data_handling(
domain_size,
periodicity=(True,) * stencil.D,
parallel=False,
default_target=target,
)
pdfs_a = dh.add_array("pdfs_a", values_per_cell=stencil.Q)
pdfs_b = dh.add_array("pdfs_b", values_per_cell=stencil.Q)
dh.fill(pdfs_a.name, 0.0, ghost_layers=True)
dh.fill(pdfs_b.name, 0.0, ghost_layers=True)
for q in range(stencil.Q):
sl = ps.make_slice[:4, :4, q] if stencil.D == 2 else ps.make_slice[:4, :4, :4, q]
dh.cpu_arrays[pdfs_a.name][sl] = q
dh.cpu_arrays[pdfs_b.name][sl] = q
dh.all_to_gpu()
direct_copy = LBMPeriodicityHandling(stencil, dh, pdfs_a.name, streaming_pattern, cupy_direct_copy=True)
kernels_copy = LBMPeriodicityHandling(stencil, dh, pdfs_b.name, streaming_pattern, cupy_direct_copy=False)
direct_copy(Timestep.EVEN)
kernels_copy(Timestep.EVEN)
dh.all_to_cpu()
np.testing.assert_equal(
dh.cpu_arrays[pdfs_a.name],
dh.cpu_arrays[pdfs_b.name]
)
@pytest.mark.parametrize(
"stencil_name", [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]
)
def test_optimised_and_full_communication_equivalence(stencil_name): def test_optimised_and_full_communication_equivalence(stencil_name):
target = ps.Target.CPU target = ps.Target.CPU
stencil = LBStencil(stencil_name) stencil = LBStencil(stencil_name)
domain_size = (4, ) * stencil.D domain_size = (4,) * stencil.D
dh = ps.create_data_handling(domain_size, periodicity=(True, ) * stencil.D, dh = ps.create_data_handling(
parallel=False, default_target=target) domain_size,
periodicity=(True,) * stencil.D,
parallel=False,
default_target=target,
)
pdf = dh.add_array("pdf", values_per_cell=len(stencil), dtype=np.int64) pdf = dh.add_array("pdf", values_per_cell=len(stencil), dtype=np.int64)
dh.fill("pdf", 0, ghost_layers=True) dh.fill("pdf", 0, ghost_layers=True)
...@@ -82,34 +195,40 @@ def test_optimised_and_full_communication_equivalence(stencil_name): ...@@ -82,34 +195,40 @@ def test_optimised_and_full_communication_equivalence(stencil_name):
gl = dh.ghost_layers_of_field("pdf") gl = dh.ghost_layers_of_field("pdf")
num = 0 num = 0
for idx, x in np.ndenumerate(dh.cpu_arrays['pdf']): for idx, x in np.ndenumerate(dh.cpu_arrays["pdf"]):
dh.cpu_arrays['pdf'][idx] = num dh.cpu_arrays["pdf"][idx] = num
dh.cpu_arrays['pdf_tmp'][idx] = num dh.cpu_arrays["pdf_tmp"][idx] = num
num += 1 num += 1
lbm_config = LBMConfig(stencil=stencil, kernel_type="stream_pull_only") lbm_config = LBMConfig(stencil=stencil, kernel_type="stream_pull_only")
lbm_opt = LBMOptimisation(symbolic_field=pdf, symbolic_temporary_field=pdf_tmp) lbm_opt = LBMOptimisation(symbolic_field=pdf, symbolic_temporary_field=pdf_tmp)
config = ps.CreateKernelConfig(target=dh.default_target, cpu_openmp=True) config = ps.CreateKernelConfig(target=dh.default_target, data_type=np.int64, cpu_openmp=True)
ac = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt) ac = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
ast = ps.create_kernel(ac, config=config) ast = ps.create_kernel(ac, config=config)
stream = ast.compile() stream = ast.compile()
full_communication = dh.synchronization_function(pdf.name, target=dh.default_target, optimization={"openmp": True}) full_communication = dh.synchronization_function(
pdf.name, target=dh.default_target, optimization={"openmp": True}
)
full_communication() full_communication()
dh.run_kernel(stream) dh.run_kernel(stream)
dh.swap("pdf", "pdf_tmp") dh.swap("pdf", "pdf_tmp")
pdf_full_communication = np.copy(dh.cpu_arrays['pdf']) pdf_full_communication = np.copy(dh.cpu_arrays["pdf"])
num = 0 num = 0
for idx, x in np.ndenumerate(dh.cpu_arrays['pdf']): for idx, x in np.ndenumerate(dh.cpu_arrays["pdf"]):
dh.cpu_arrays['pdf'][idx] = num dh.cpu_arrays["pdf"][idx] = num
dh.cpu_arrays['pdf_tmp'][idx] = num dh.cpu_arrays["pdf_tmp"][idx] = num
num += 1 num += 1
optimised_communication = LBMPeriodicityHandling(stencil=stencil, data_handling=dh, pdf_field_name=pdf.name, optimised_communication = LBMPeriodicityHandling(
streaming_pattern='pull') stencil=stencil,
data_handling=dh,
pdf_field_name=pdf.name,
streaming_pattern="pull",
)
optimised_communication() optimised_communication()
dh.run_kernel(stream) dh.run_kernel(stream)
dh.swap("pdf", "pdf_tmp") dh.swap("pdf", "pdf_tmp")
...@@ -119,9 +238,14 @@ def test_optimised_and_full_communication_equivalence(stencil_name): ...@@ -119,9 +238,14 @@ def test_optimised_and_full_communication_equivalence(stencil_name):
for j in range(gl, domain_size[1]): for j in range(gl, domain_size[1]):
for k in range(gl, domain_size[2]): for k in range(gl, domain_size[2]):
for f in range(len(stencil)): for f in range(len(stencil)):
assert dh.cpu_arrays['pdf'][i, j, k, f] == pdf_full_communication[i, j, k, f], print(f) assert (
dh.cpu_arrays["pdf"][i, j, k, f]
== pdf_full_communication[i, j, k, f]
), print(f)
else: else:
for i in range(gl, domain_size[0]): for i in range(gl, domain_size[0]):
for j in range(gl, domain_size[1]): for j in range(gl, domain_size[1]):
for f in range(len(stencil)): for f in range(len(stencil)):
assert dh.cpu_arrays['pdf'][i, j, f] == pdf_full_communication[i, j, f] assert (
dh.cpu_arrays["pdf"][i, j, f] == pdf_full_communication[i, j, f]
)
...@@ -2,9 +2,16 @@ from statistics import median ...@@ -2,9 +2,16 @@ from statistics import median
import numpy as np import numpy as np
import sympy as sp import sympy as sp
import pytest
from lbmpy.scenarios import create_lid_driven_cavity from lbmpy.scenarios import create_lid_driven_cavity
from pystencils.cpu.cpujit import add_or_change_compiler_flags from lbmpy._compat import IS_PYSTENCILS_2
if IS_PYSTENCILS_2:
pytest.skip(reason="Benchmark test case needs to be rebuilt for pystencils 2.0", allow_module_level=True)
else:
from pystencils.cpu.cpujit import add_or_change_compiler_flags
from pystencils.runhelper import ParameterStudy from pystencils.runhelper import ParameterStudy
......
...@@ -58,7 +58,12 @@ def flow_around_sphere(stencil, galilean_correction, fourth_order_correction, L_ ...@@ -58,7 +58,12 @@ def flow_around_sphere(stencil, galilean_correction, fourth_order_correction, L_
boundary_assignments) boundary_assignments)
iter_slice = get_slice_before_ghost_layer((1,) + (0,) * (stencil.D - 1)) iter_slice = get_slice_before_ghost_layer((1,) + (0,) * (stencil.D - 1))
extrapolation_ast = create_kernel( extrapolation_ast = create_kernel(
boundary_assignments, config=CreateKernelConfig(iteration_slice=iter_slice, ghost_layers=1, target=target)) boundary_assignments,
config=CreateKernelConfig(
iteration_slice=iter_slice,
target=target,
skip_independence_check=True)
)
return extrapolation_ast.compile() return extrapolation_ast.compile()
dh = create_data_handling(channel_size, periodicity=False, default_layout='fzyx', default_target=target) dh = create_data_handling(channel_size, periodicity=False, default_layout='fzyx', default_target=target)
......
...@@ -10,7 +10,9 @@ import warnings ...@@ -10,7 +10,9 @@ import warnings
import numpy as np import numpy as np
import pytest import pytest
from pystencils.backends.simd_instruction_sets import get_supported_instruction_sets from pystencils import Target
from lbmpy._compat import get_supported_instruction_sets
from lbmpy.boundaries.boundaryconditions import NoSlip from lbmpy.boundaries.boundaryconditions import NoSlip
from lbmpy.geometry import get_pipe_velocity_field from lbmpy.geometry import get_pipe_velocity_field
from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
......
...@@ -202,19 +202,19 @@ def create_full_parameter_study(): ...@@ -202,19 +202,19 @@ def create_full_parameter_study():
config = CreateKernelConfig(target=Target.CPU) config = CreateKernelConfig(target=Target.CPU)
methods = [LBMConfig(method=Method.TRT, relaxation_rates=[omega, 1]), methods = [LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT, relaxation_rates=[omega, 1]),
LBMConfig(method=Method.MRT, relaxation_rates=[omega], equilibrium_order=2), LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.MRT, relaxation_rates=[omega],
LBMConfig(method=Method.MRT, relaxation_rates=[omega], equilibrium_order=3), equilibrium_order=2),
LBMConfig(method=Method.CUMULANT, relaxation_rates=[omega], # cumulant LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.MRT, relaxation_rates=[omega],
equilibrium_order=3),
LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, relaxation_rates=[omega], # cumulant
compressible=True, continuous_equilibrium=True, equilibrium_order=3), compressible=True, continuous_equilibrium=True, equilibrium_order=3),
LBMConfig(method=Method.CUMULANT, relaxation_rates=[omega], # cumulant with fourth-order correction LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, relaxation_rates=[omega], # cumulant with fourth-order correction
compressible=True, continuous_equilibrium=True, fourth_order_correction=0.1), compressible=True, continuous_equilibrium=True, fourth_order_correction=0.1),
LBMConfig(method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f], entropic=True, LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f],
zero_centered=False, # entropic order 2 entropic=True, zero_centered=False, continuous_equilibrium=True, equilibrium_order=2),
continuous_equilibrium=True, equilibrium_order=2), LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f],
LBMConfig(method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f], entropic=True, entropic=True, zero_centered=False, continuous_equilibrium=True, equilibrium_order=3),
zero_centered=False, # entropic order 3
continuous_equilibrium=True, equilibrium_order=3),
# entropic cumulant: not supported for the moment # entropic cumulant: not supported for the moment
# LBMConfig(method=Method.CUMULANT, relaxation_rates=["omega", "omega_f"], entropic=True, zero_centered=False, # LBMConfig(method=Method.CUMULANT, relaxation_rates=["omega", "omega_f"], entropic=True, zero_centered=False,
......
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import pytest import pytest
pytest.importorskip('cupy') pytest.importorskip('cupy')
``` ```
%% Output %% Output
<module 'cupy' from '/home/markus/.local/lib/python3.11/site-packages/cupy/__init__.py'> <module 'cupy' from '/home/markus/.local/lib/python3.11/site-packages/cupy/__init__.py'>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from lbmpy.session import * from lbmpy.session import *
from lbmpy.phasefield.n_phase_boyer import * from lbmpy.phasefield.n_phase_boyer import *
from lbmpy.phasefield.kerneleqs import * from lbmpy.phasefield.kerneleqs import *
from lbmpy.phasefield.contact_angle_circle_fitting import * from lbmpy.phasefield.contact_angle_circle_fitting import *
from scipy.ndimage.filters import gaussian_filter from scipy.ndimage import gaussian_filter
from pystencils.simp import sympy_cse_on_assignment_list from pystencils.simp import sympy_cse_on_assignment_list
one = sp.sympify(1) one = sp.sympify(1)
import pyximport import pyximport
pyximport.install(language_level=3) pyximport.install(language_level=3)
from lbmpy.phasefield.simplex_projection import simplex_projection_2d # NOQA from lbmpy.phasefield.simplex_projection import simplex_projection_2d # NOQA
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Simulation arbitrary surface tension case # Simulation arbitrary surface tension case
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
n = 4 n = 4
dx, dt = 1, 1 dx, dt = 1, 1
mobility = 2e-3 mobility = 2e-3
domain_size = (150, 150) domain_size = (150, 150)
ε = one * 4 ε = one * 4
penalty_factor = 0 penalty_factor = 0
stabilization_factor = 10 stabilization_factor = 10
κ = (one, one/2, one/3, one/4) κ = (one, one/2, one/3, one/4)
sigma_factor = one / 15 sigma_factor = one / 15
σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 ) σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
dh = create_data_handling(domain_size, periodicity=True, default_target=ps.Target.GPU) dh = create_data_handling(domain_size, periodicity=True, default_target=ps.Target.GPU)
c = dh.add_array('c', values_per_cell=n) c = dh.add_array('c', values_per_cell=n)
c_tmp = dh.add_array_like('c_tmp', 'c') c_tmp = dh.add_array_like('c_tmp', 'c')
μ = dh.add_array('mu', values_per_cell=n) μ = dh.add_array('mu', values_per_cell=n)
cvec = c.center_vector cvec = c.center_vector
μvec = μ.center_vector μvec = μ.center_vector
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
α, _ = diffusion_coefficients(σ) α, _ = diffusion_coefficients(σ)
f = lambda c: c**2 * ( 1 - c ) **2 f = lambda c: c**2 * ( 1 - c ) **2
a, b = compute_ab(f) a, b = compute_ab(f)
capital_f = capital_f0(cvec, σ) + correction_g(cvec, σ) + stabilization_factor * stabilization_term(cvec, α) capital_f = capital_f0(cvec, σ) + correction_g(cvec, σ) + stabilization_factor * stabilization_term(cvec, α)
f_bulk = free_energy_bulk(capital_f, b, ε) + penalty_factor * (one - sum(cvec)) f_bulk = free_energy_bulk(capital_f, b, ε) + penalty_factor * (one - sum(cvec))
f_if = free_energy_interfacial(cvec, σ, a, ε) f_if = free_energy_interfacial(cvec, σ, a, ε)
f = f_bulk + f_if f = f_bulk + f_if
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#f_bulk #f_bulk
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
μ_assignments = mu_kernel(f, cvec, c, μ) μ_assignments = mu_kernel(f, cvec, c, μ)
μ_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in μ_assignments] μ_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in μ_assignments]
μ_assignments = sympy_cse_on_assignment_list(μ_assignments) μ_assignments = sympy_cse_on_assignment_list(μ_assignments)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
discretize = fd.Discretization2ndOrder(dx=dx, dt=dt) discretize = fd.Discretization2ndOrder(dx=dx, dt=dt)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def lapl(e): def lapl(e):
return sum(ps.fd.diff(e, d, d) for d in range(dh.dim)) return sum(ps.fd.diff(e, d, d) for d in range(dh.dim))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
rhs = α * μvec rhs = α * μvec
discretized_rhs = [discretize(fd.expand_diff_full( lapl(mobility * rhs_i) + fd.transient(cvec[i], idx=i), functions=μvec)) discretized_rhs = [discretize(fd.expand_diff_full( lapl(mobility * rhs_i) + fd.transient(cvec[i], idx=i), functions=μvec))
for i, rhs_i in enumerate(rhs)] for i, rhs_i in enumerate(rhs)]
c_assignments = [Assignment(lhs, rhs) c_assignments = [Assignment(lhs, rhs)
for lhs, rhs in zip(c_tmp.center_vector, discretized_rhs)] for lhs, rhs in zip(c_tmp.center_vector, discretized_rhs)]
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#c_assignments #c_assignments
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
μ_sync = dh.synchronization_function(μ.name) μ_sync = dh.synchronization_function(μ.name)
c_sync = dh.synchronization_function(c.name) c_sync = dh.synchronization_function(c.name)
optimization = {'cpu_openmp': False, 'cpu_vectorize_info': None} optimization = {'cpu_openmp': False, 'cpu_vectorize_info': None}
config = ps.CreateKernelConfig(cpu_openmp=False, target=dh.default_target) config = ps.CreateKernelConfig(cpu_openmp=False, target=dh.default_target)
μ_kernel = create_kernel(μ_assignments, config=config).compile() μ_kernel = create_kernel(μ_assignments, config=config).compile()
c_kernel = create_kernel(c_assignments, config=config).compile() c_kernel = create_kernel(c_assignments, config=config).compile()
def set_c(slice_obj, values): def set_c(slice_obj, values):
for block in dh.iterate(slice_obj): for block in dh.iterate(slice_obj):
arr = block[c.name] arr = block[c.name]
arr[..., : ] = values arr[..., : ] = values
def smooth(): def smooth():
for block in dh.iterate(ghost_layers=True): for block in dh.iterate(ghost_layers=True):
c_arr = block[c.name] c_arr = block[c.name]
for i in range(n): for i in range(n):
gaussian_filter(c_arr[..., i], sigma=2, output=c_arr[..., i]) gaussian_filter(c_arr[..., i], sigma=2, output=c_arr[..., i])
def time_loop(steps): def time_loop(steps):
dh.all_to_gpu() dh.all_to_gpu()
for t in range(steps): for t in range(steps):
c_sync() c_sync()
dh.run_kernel(μ_kernel) dh.run_kernel(μ_kernel)
μ_sync() μ_sync()
dh.run_kernel(c_kernel) dh.run_kernel(c_kernel)
dh.swap(c.name, c_tmp.name) dh.swap(c.name, c_tmp.name)
#simplex_projection_2d(dh.cpu_arrays[c.name]) #simplex_projection_2d(dh.cpu_arrays[c.name])
dh.all_to_cpu() dh.all_to_cpu()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
set_c(make_slice[:, :], [0, 0, 0, 0]) set_c(make_slice[:, :], [0, 0, 0, 0])
set_c(make_slice[:, 0.5:], [1, 0, 0, 0]) set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
set_c(make_slice[:, :0.5], [0, 1, 0, 0]) set_c(make_slice[:, :0.5], [0, 1, 0, 0])
set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0]) set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
smooth() smooth()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#dh.load_all('n_phases_state_size200_stab10.npz') #dh.load_all('n_phases_state_size200_stab10.npz')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
plt.phase_plot(dh.gather_array(c.name)) plt.phase_plot(dh.gather_array(c.name))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j])) neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j]))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import time import time
for i in range(10): for i in range(10):
start = time.perf_counter() start = time.perf_counter()
time_loop(1_000) time_loop(1_000)
end = time.perf_counter() end = time.perf_counter()
try: try:
print(i, end - start, liquid_lens_neumann_angles(dh.gather_array(c.name))) print(i, end - start, liquid_lens_neumann_angles(dh.gather_array(c.name)))
except Exception: except Exception:
print(i, end - start, "none found") print(i, end - start, "none found")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
plt.subplot(1,3,1) plt.subplot(1,3,1)
t = dh.gather_array(c.name, make_slice[25, :]).squeeze() t = dh.gather_array(c.name, make_slice[25, :]).squeeze()
plt.plot(t); plt.plot(t);
plt.subplot(1,3,2) plt.subplot(1,3,2)
plt.phase_plot(dh.gather_array(c.name), linewidth=1) plt.phase_plot(dh.gather_array(c.name), linewidth=1)
plt.subplot(1,3,3) plt.subplot(1,3,3)
plt.scalar_field(dh.gather_array(μ.name)[:, :, 2]) plt.scalar_field(dh.gather_array(μ.name)[:, :, 2])
plt.colorbar(); plt.colorbar();
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
assert not np.isnan(dh.max(c.name)) assert not np.isnan(dh.max(c.name))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
t = dh.gather_array(c.name, make_slice[25, 55:90]).squeeze() t = dh.gather_array(c.name, make_slice[25, 55:90]).squeeze()
plt.hlines(0.5, 0, 30) plt.hlines(0.5, 0, 30)
plt.plot(t); plt.plot(t);
``` ```
......
...@@ -44,6 +44,6 @@ def test_analytical(): ...@@ -44,6 +44,6 @@ def test_analytical():
assert np.isclose(parameters.density_heavy, 1.0) assert np.isclose(parameters.density_heavy, 1.0)
assert np.isclose(parameters.density_light, 0.001207114228456914) assert np.isclose(parameters.density_light, 0.001207114228456914)
assert np.isclose(parameters.dynamic_viscosity_heavy, 5.733727652152216e-05) assert np.isclose(parameters.dynamic_viscosity_heavy, 5.733727652152216e-05)
assert np.isclose(parameters.dynamic_viscosity_light, 0.0008630017037694861) assert np.isclose(parameters.dynamic_viscosity_light, 1.0417416358027054e-06)
assert np.isclose(parameters.gravitational_acceleration, -7.407407407407407e-08) assert np.isclose(parameters.gravitational_acceleration, -7.407407407407407e-08)
assert np.isclose(parameters.surface_tension, 3.149857262258028e-05, rtol=1e-05) assert np.isclose(parameters.surface_tension, 3.149857262258028e-05, rtol=1e-05)
import math import math
import pytest
import pystencils as ps import pystencils as ps
from pystencils.boundaries.boundaryhandling import BoundaryHandling from pystencils.boundaries.boundaryhandling import BoundaryHandling
from lbmpy.enums import Stencil from lbmpy.enums import Stencil
from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
from lbmpy.stencils import LBStencil from lbmpy.stencils import LBStencil
from lbmpy._compat import IS_PYSTENCILS_2
import numpy as np import numpy as np
@pytest.mark.skipif(
IS_PYSTENCILS_2,
reason="Contact angle calculation not yet available with pystencils 2.0"
)
def test_contact_angle(): def test_contact_angle():
from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
stencil = LBStencil(Stencil.D2Q9) stencil = LBStencil(Stencil.D2Q9)
contact_angle = 45 contact_angle = 45
phase_value = 0.5 phase_value = 0.5
......
import numpy as np import numpy as np
import pytest import pytest
from lbmpy.boundaries import NoSlip, UBB, SimpleExtrapolationOutflow, ExtrapolationOutflow, \ from lbmpy.boundaries import (NoSlip, NoSlipLinearBouzidi, QuadraticBounceBack,
FixedDensity, DiffusionDirichlet, NeumannByCopy, StreamInConstant, FreeSlip UBB, SimpleExtrapolationOutflow, ExtrapolationOutflow, FixedDensity, DiffusionDirichlet,
from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling NeumannByCopy, StreamInConstant, FreeSlip)
from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling, create_lattice_boltzmann_boundary_kernel
from lbmpy.creationfunctions import create_lb_function, create_lb_method, LBMConfig from lbmpy.creationfunctions import create_lb_function, create_lb_method, LBMConfig
from lbmpy.enums import Stencil, Method from lbmpy.enums import Stencil, Method
from lbmpy.geometry import add_box_boundary from lbmpy.geometry import add_box_boundary
from lbmpy.lbstep import LatticeBoltzmannStep from lbmpy.lbstep import LatticeBoltzmannStep
from lbmpy.stencils import LBStencil from lbmpy.stencils import LBStencil
import pystencils as ps
from pystencils import create_data_handling, make_slice, Target, CreateKernelConfig from pystencils import create_data_handling, make_slice, Target, CreateKernelConfig
from pystencils.slicing import slice_from_direction from pystencils.slicing import slice_from_direction
from pystencils.stencil import inverse_direction from pystencils.stencil import inverse_direction
from lbmpy._compat import IS_PYSTENCILS_2
def mirror_stencil(direction, mirror_axis): def mirror_stencil(direction, mirror_axis):
for i, n in enumerate(mirror_axis): for i, n in enumerate(mirror_axis):
...@@ -435,3 +439,61 @@ def test_boundary_utility_functions(): ...@@ -435,3 +439,61 @@ def test_boundary_utility_functions():
assert stream == StreamInConstant(constant=1.0, name="stream") assert stream == StreamInConstant(constant=1.0, name="stream")
assert not stream == StreamInConstant(constant=1.0, name="test") assert not stream == StreamInConstant(constant=1.0, name="test")
assert not stream == noslip assert not stream == noslip
@pytest.mark.parametrize("given_force_vector", [True, False])
@pytest.mark.parametrize("dtype", ["float32", "float64"])
def test_force_on_boundary(given_force_vector, dtype):
stencil = LBStencil(Stencil.D2Q9)
pdfs = ps.fields(f"pdfs_src({stencil.Q}): {dtype}[{stencil.D}D]", layout='fzyx')
method = create_lb_method(lbm_config=LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=1.8))
noslip = NoSlip(name="noslip", calculate_force_on_boundary=True)
bouzidi = NoSlipLinearBouzidi(
name="bouzidi",
data_type=dtype,
calculate_force_on_boundary=True
)
qq_bounce_Back = QuadraticBounceBack(
name="qqBB",
relaxation_rate=1.8,
data_type=dtype,
calculate_force_on_boundary=True
)
boundary_objects = [noslip, bouzidi, qq_bounce_Back]
for boundary in boundary_objects:
if given_force_vector:
force_vector_type = np.dtype([(f"F_{i}", dtype) for i in range(stencil.D)], align=True)
force_vector = ps.Field('forceVector', ps.FieldType.INDEXED, force_vector_type, layout=[0],
shape=(ps.TypedSymbol("forceVectorSize", "int32"), 1), strides=(1, 1))
else:
force_vector = None
index_struct_dtype = _numpy_data_type_for_boundary_object(boundary, stencil.D)
index_field = ps.Field('indexVector', ps.FieldType.INDEXED, index_struct_dtype, layout=[0],
shape=(ps.TypedSymbol("indexVectorSize", "int32"), 1), strides=(1, 1))
create_lattice_boltzmann_boundary_kernel(
pdfs,
index_field,
method,
boundary,
force_vector=force_vector,
**({"default_dtype": dtype} if IS_PYSTENCILS_2 else dict())
)
def _numpy_data_type_for_boundary_object(boundary_object, dim):
boundary_index_array_coordinate_names = ["x", "y", "z"]
direction_member_name = "dir"
default_index_array_dtype = np.int32
coordinate_names = boundary_index_array_coordinate_names[:dim]
return np.dtype(
[(name, default_index_array_dtype) for name in coordinate_names]
+ [(direction_member_name, default_index_array_dtype)]
+ [(i[0], i[1].numpy_dtype) for i in boundary_object.additional_data],
align=True,)
This diff is collapsed.
...@@ -2,19 +2,23 @@ from dataclasses import replace ...@@ -2,19 +2,23 @@ from dataclasses import replace
import numpy as np import numpy as np
import pytest import pytest
from pystencils import Backend, Target, CreateKernelConfig from pystencils import Target, CreateKernelConfig
from lbmpy.creationfunctions import LBMConfig, LBMOptimisation from lbmpy.creationfunctions import LBMConfig, LBMOptimisation
from lbmpy.enums import Method, Stencil from lbmpy.enums import Method, Stencil
from lbmpy.scenarios import create_channel from lbmpy.scenarios import create_channel
from lbmpy.stencils import LBStencil from lbmpy.stencils import LBStencil
from lbmpy._compat import IS_PYSTENCILS_2
def run_equivalence_test(domain_size, lbm_config, lbm_opt, config, time_steps=13): def run_equivalence_test(domain_size, lbm_config, lbm_opt, base_config, time_steps=13):
config = replace(config, target=Target.CPU) config = replace(base_config, target=Target.CPU)
cpu_scenario = create_channel(domain_size=domain_size, pressure_difference=0.001, cpu_scenario = create_channel(domain_size=domain_size, pressure_difference=0.001,
lbm_config=lbm_config, lbm_optimisation=lbm_opt, config=config) lbm_config=lbm_config, lbm_optimisation=lbm_opt, config=config)
config = replace(config, target=Target.GPU, backend=Backend.CUDA) config = replace(base_config, target=Target.GPU)
if not IS_PYSTENCILS_2:
from pystencils.enums import Backend
config = replace(config, backend=Backend.CUDA)
gpu_scenario = create_channel(domain_size=domain_size, pressure_difference=0.001, gpu_scenario = create_channel(domain_size=domain_size, pressure_difference=0.001,
lbm_config=lbm_config, lbm_optimisation=lbm_opt, config=config) lbm_config=lbm_config, lbm_optimisation=lbm_opt, config=config)
...@@ -48,6 +52,10 @@ def test_force_driven_channel_short(scenario): ...@@ -48,6 +52,10 @@ def test_force_driven_channel_short(scenario):
if block_size is not False: if block_size is not False:
config = CreateKernelConfig(gpu_indexing_params={'block_size': block_size}) config = CreateKernelConfig(gpu_indexing_params={'block_size': block_size})
else: else:
config = CreateKernelConfig(gpu_indexing='line') if IS_PYSTENCILS_2:
config = CreateKernelConfig()
config.gpu.indexing_scheme = "blockwise4d"
else:
config = CreateKernelConfig(gpu_indexing='line')
run_equivalence_test(domain_size=ds, lbm_config=lbm_config, lbm_opt=lbm_opt, config=config) run_equivalence_test(domain_size=ds, lbm_config=lbm_config, lbm_opt=lbm_opt, base_config=config)