Commit 6cf0ce87 authored by Markus Holzer's avatar Markus Holzer
Browse files

Implemented short test case

parent 730c62d4
Pipeline #32562 passed with stage
in 32 minutes and 11 seconds
......@@ -126,7 +126,6 @@ class UBB(LbBoundary):
"""
def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None, data_type='double'):
super(UBB, self).__init__(name)
self._velocity = velocity
self._adaptVelocityToForce = adapt_velocity_to_force
if callable(self._velocity) and not dim:
......@@ -136,12 +135,14 @@ class UBB(LbBoundary):
self.dim = dim
self.data_type = data_type
super(UBB, self).__init__(name)
@property
def additional_data(self):
""" In case of the UBB boundary additional data is a velocity vector. This vector is added to each cell to
realize velocity profiles for the inlet."""
if self.velocity_is_callable:
return [('vel_%d' % (i,), create_type(self.data_type)) for i in range(self.dim)]
return [(f'vel_{i}', create_type(self.data_type)) for i in range(self.dim)]
else:
return []
......@@ -211,8 +212,6 @@ class UBB(LbBoundary):
else:
return [Assignment(f_in(inv_dir[direction]),
f_out(direction) - vel_term)]
# end class UBB
......@@ -294,11 +293,10 @@ class ExtrapolationOutflow(LbBoundary):
streaming_pattern='pull', zeroth_timestep=Timestep.BOTH,
initial_density=None, initial_velocity=None, data_type='double'):
self.data_type = data_type
self.lb_method = lb_method
self.stencil = lb_method.stencil
self.dim = len(self.stencil[0])
if isinstance(normal_direction, str):
normal_direction = direction_string_to_offset(normal_direction, dim=self.dim)
......@@ -316,6 +314,8 @@ class ExtrapolationOutflow(LbBoundary):
self.initial_velocity = initial_velocity
self.equilibrium_calculation = None
self.data_type = data_type
if initial_density and initial_velocity:
equilibrium = lb_method.get_equilibrium(conserved_quantity_equations=AssignmentCollection([]))
rho = lb_method.zeroth_order_equilibrium_moment_symbol
......@@ -405,7 +405,6 @@ class ExtrapolationOutflow(LbBoundary):
return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
# end class ExtrapolationOutflow
......@@ -420,8 +419,9 @@ class FixedDensity(LbBoundary):
def __init__(self, density, name=None):
if name is None:
name = "Fixed Density " + str(density)
self.density = density
super(FixedDensity, self).__init__(name)
self._density = density
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
def remove_asymmetric_part_of_main_assignments(assignment_collection, degrees_of_freedom):
......@@ -441,7 +441,7 @@ class FixedDensity(LbBoundary):
density_symbol = cqc.defined_symbols()['density']
density = self._density
density = self.density
equilibrium_input = cqc.equilibrium_input_equations_from_init_values(density=density)
equilibrium_input = equilibrium_input.new_without_subexpressions()
density_eq = equilibrium_input.main_assignments[0]
......@@ -473,8 +473,9 @@ class DiffusionDirichlet(LbBoundary):
def __init__(self, concentration, name=None):
if name is None:
name = "Diffusion Dirichlet " + str(concentration)
self.concentration = concentration
super(DiffusionDirichlet, self).__init__(name)
self._concentration = concentration
def get_additional_code_nodes(self, lb_method):
"""Return a list of code nodes that will be added in the generated code before the index field loop.
......@@ -490,7 +491,7 @@ class DiffusionDirichlet(LbBoundary):
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
w_dir = LbmWeightInfo.weight_of_direction(dir_symbol, lb_method)
return [Assignment(f_in(inv_dir[dir_symbol]),
2 * w_dir * self._concentration - f_out(dir_symbol))]
2 * w_dir * self.concentration - f_out(dir_symbol))]
# end class DiffusionDirichlet
......@@ -520,10 +521,6 @@ class NeumannByCopy(LbBoundary):
# All boundaries of these class behave equal -> should also be equal
return hash("NeumannByCopy")
def __eq__(self, other):
return type(other) == NeumannByCopy
# end class NeumannByCopy
......@@ -538,7 +535,7 @@ class StreamInConstant(LbBoundary):
def __init__(self, constant, name=None):
super(StreamInConstant, self).__init__(name)
self._constant = constant
self.constant = constant
def get_additional_code_nodes(self, lb_method):
"""Return a list of code nodes that will be added in the generated code before the index field loop.
......@@ -553,13 +550,11 @@ class StreamInConstant(LbBoundary):
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
return [Assignment(f_in(inv_dir[dir_symbol]), self._constant),
Assignment(f_out[neighbour_offset](dir_symbol), self._constant)]
return [Assignment(f_in(inv_dir[dir_symbol]), self.constant),
Assignment(f_out[neighbour_offset](dir_symbol), self.constant)]
def __hash__(self):
# All boundaries of these class behave equal -> should also be equal
return hash("StreamInConstant")
def __eq__(self, other):
return type(other) == StreamInConstant
# end class StreamInConstant
......@@ -12,7 +12,7 @@ from sympy import Rational as R
from pystencils.cache import disk_cache
def get_weights(stencil, c_s_sq):
def get_weights(stencil, c_s_sq=sp.Rational(1, 3)):
q = len(stencil)
if c_s_sq != sp.Rational(1, 3) and c_s_sq != sp.Symbol("c_s") ** 2:
......
from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
from lbmpy.session import *
from lbmpy.stencils import get_stencil
import math
import pytest
def test_diffusion():
"""
Runs the "Diffusion from Plate in Uniform Flow" benchmark as it is described
in [ch. 8.6.3, The Lattice Boltzmann Method, Krüger et al.].
import pystencils as ps
from pystencils.slicing import slice_from_direction
from lbmpy import pdf_initialization_assignments
from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
from lbmpy.boundaries.boundaryconditions import DiffusionDirichlet, NeumannByCopy
from lbmpy.geometry import add_box_boundary
from lbmpy.stencils import get_stencil
from lbmpy.maxwellian_equilibrium import get_weights
import sympy as sp
import numpy as np
def test_diffusion_boundary():
domain_size = (10, 10)
stencil = get_stencil('D2Q9')
weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
concentration = 1.0
# Data Handling
dh = ps.create_data_handling(domain_size=domain_size)
dC/dy = 0
┌───────────────┐
│ → → → │
C = 0 │ → u → │ dC/dx = 0
│ → → → │
└───────────────┘
C = 1
dh.add_array('pdfs', values_per_cell=len(stencil))
dh.fill("pdfs", 0.0, ghost_layers=True)
The analytical solution is given by:
C(x,y) = 1 * erfc(y / sqrt(4Dx/u))
method = create_lb_method(stencil=stencil, method="srt", relaxation_rate=1.8, compressible=True)
The hydrodynamic field is not simulated, instead a constant velocity is assumed.
# Boundary Handling
bh = LatticeBoltzmannBoundaryHandling(method, dh, 'pdfs', name="bh")
add_box_boundary(bh, boundary=DiffusionDirichlet(concentration))
bh()
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[0, 1:-2, 4] - 2 * weights[4]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[0, 1:-2, 6] - 2 * weights[6]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[0, 2:-1, 8] - 2 * weights[8]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[-1, 1:-2, 3] - 2 * weights[3]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[-1, 1:-2, 5] - 2 * weights[5]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[-1, 2:-1, 7] - 2 * weights[7]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[1:-2, 0, 1] - 2 * weights[1]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[2:, 0, 5] - 2 * weights[5]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[:-2, 0, 6] - 2 * weights[6]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[1:-2, -1, 2] - 2 * weights[2]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[2:, -1, 7] - 2 * weights[7]) < 1e-14)
assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[:-2, -1, 8] - 2 * weights[8]) < 1e-14)
@pytest.mark.longrun
def test_diffusion():
"""
Runs the "Diffusion from Plate in Uniform Flow" benchmark as it is described
in [ch. 8.6.3, The Lattice Boltzmann Method, Krüger et al.].
dC/dy = 0
┌───────────────┐
│ → → → │
C = 0 │ → u → │ dC/dx = 0
│ → → → │
└───────────────┘
C = 1
The analytical solution is given by:
C(x,y) = 1 * erfc(y / sqrt(4Dx/u))
The hydrodynamic field is not simulated, instead a constant velocity is assumed.
"""
# Parameters
domain_size = (1600, 160)
omega = 1.38
D = (1/omega - 0.5) / 3
velocity = 0.05
time_steps = 50000
stencil = get_stencil('D2Q9')
# Data Handling
dh = ps.create_data_handling(domain_size=domain_size)
vel_field = dh.add_array('vel_field', values_per_cell=dh.dim)
dh.fill('vel_field', velocity, 0, ghost_layers=True)
dh.fill('vel_field', 0.0 , 1, ghost_layers=True)
con_field = dh.add_array('con_field', values_per_cell=1)
dh.fill('con_field', 0.0, ghost_layers=True)
pdfs = dh.add_array('pdfs' , values_per_cell=len(stencil))
pdfs_tmp = dh.add_array('pdfs_tmp', values_per_cell=len(stencil))
# Lattice Boltzmann method
params = { 'stencil' : stencil,
'method' : 'mrt',
'relaxation_rates' : [1, 1.5, 1, 1.5, 1],
'velocity_input' : vel_field,
'output' : {'density' : con_field},
'compressible' : True,
'weighted' : True,
'optimization' : {'symbolic_field' : pdfs,
'symbolic_temporary_field' : pdfs_tmp},
'kernel_type' : 'stream_pull_collide'
}
method = create_lb_method(**params)
method.set_conserved_moments_relaxation_rate(omega)
update_rule = create_lb_update_rule(lb_method=method, **params)
kernel = ps.create_kernel(update_rule, target=dh.default_target).compile()
# PDF initalization
init = pdf_initialization_assignments(method, con_field.center,
vel_field.center_vector, pdfs.center_vector)
dh.run_kernel(ps.create_kernel(init, target=dh.default_target).compile())
# Boundary Handling
bh = LatticeBoltzmannBoundaryHandling(update_rule.method, dh, 'pdfs', name="bh")
add_box_boundary(bh, boundary=NeumannByCopy())
bh.set_boundary(DiffusionDirichlet(0), slice_from_direction('W', dh.dim))
bh.set_boundary(DiffusionDirichlet(1), slice_from_direction('S', dh.dim))
# Timeloop
for i in range(time_steps):
bh()
dh.run_kernel(kernel)
dh.swap("pdfs", "pdfs_tmp")
# Verification
x = np.arange(1, domain_size[0], 1)
y = np.arange(0, domain_size[1], 1)
X, Y = np.meshgrid(x, y)
analytical = np.zeros(domain_size)
analytical[1:,:] = np.vectorize(math.erfc)(Y / np.vectorize(math.sqrt)(4 * D * X / velocity)).transpose()
simulated = dh.gather_array('con_field', ghost_layers=False)
residual = 0
for i in x:
for j in y:
residual += (simulated[i,j] - analytical[i,j])**2
residual = math.sqrt(residual / (domain_size[0] * domain_size[1]))
assert residual < 1e-2
# Parameters
domain_size = (1600, 160)
omega = 1.38
diffusion = (1 / omega - 0.5) / 3
velocity = 0.05
time_steps = 50000
stencil = get_stencil('D2Q9')
# Data Handling
dh = ps.create_data_handling(domain_size=domain_size)
vel_field = dh.add_array('vel_field', values_per_cell=dh.dim)
dh.fill('vel_field', velocity, 0, ghost_layers=True)
dh.fill('vel_field', 0.0, 1, ghost_layers=True)
con_field = dh.add_array('con_field', values_per_cell=1)
dh.fill('con_field', 0.0, ghost_layers=True)
pdfs = dh.add_array('pdfs', values_per_cell=len(stencil))
pdfs_tmp = dh.add_array('pdfs_tmp', values_per_cell=len(stencil))
# Lattice Boltzmann method
params = {'stencil': stencil,
'method': 'mrt',
'relaxation_rates': [1, 1.5, 1, 1.5, 1],
'velocity_input': vel_field,
'output': {'density': con_field},
'compressible': True,
'weighted': True,
'optimization': {'symbolic_field': pdfs,
'symbolic_temporary_field': pdfs_tmp},
'kernel_type': 'stream_pull_collide'
}
method = create_lb_method(**params)
method.set_conserved_moments_relaxation_rate(omega)
update_rule = create_lb_update_rule(lb_method=method, **params)
kernel = ps.create_kernel(update_rule).compile()
# PDF initalization
init = pdf_initialization_assignments(method, con_field.center,
vel_field.center_vector, pdfs.center_vector)
dh.run_kernel(ps.create_kernel(init).compile())
# Boundary Handling
bh = LatticeBoltzmannBoundaryHandling(update_rule.method, dh, 'pdfs', name="bh")
add_box_boundary(bh, boundary=NeumannByCopy())
bh.set_boundary(DiffusionDirichlet(0), slice_from_direction('W', dh.dim))
bh.set_boundary(DiffusionDirichlet(1), slice_from_direction('S', dh.dim))
# Timeloop
for i in range(time_steps):
bh()
dh.run_kernel(kernel)
dh.swap("pdfs", "pdfs_tmp")
# Verification
x = np.arange(1, domain_size[0], 1)
y = np.arange(0, domain_size[1], 1)
X, Y = np.meshgrid(x, y)
analytical = np.zeros(domain_size)
analytical[1:, :] = np.vectorize(math.erfc)(Y / np.vectorize(math.sqrt)(4 * diffusion * X / velocity)).transpose()
simulated = dh.gather_array('con_field', ghost_layers=False)
residual = 0
for i in x:
for j in y:
residual += (simulated[i, j] - analytical[i, j]) ** 2
residual = math.sqrt(residual / (domain_size[0] * domain_size[1]))
assert residual < 1e-2
......@@ -15,12 +15,13 @@ def test_lbm_vectorization_short():
ldc1_ref.run(10)
ldc1 = create_lid_driven_cavity(size1, relaxation_rate=relaxation_rate,
optimization={'vectorization': {'instruction_set': get_supported_instruction_sets()[-1],
'assume_aligned': True,
'nontemporal': True,
'assume_inner_stride_one': True,
'assume_sufficient_line_padding': False,
}},
optimization={
'vectorization': {'instruction_set': get_supported_instruction_sets()[-1],
'assume_aligned': True,
'nontemporal': True,
'assume_inner_stride_one': True,
'assume_sufficient_line_padding': False,
}},
fixed_loop_sizes=False)
ldc1.run(10)
......@@ -33,10 +34,10 @@ def test_lbm_vectorization_short():
@pytest.mark.longrun
def test_lbm_vectorization(instruction_set, aligned_and_padding, nontemporal, double_precision, fixed_loop_sizes):
vectorization_options = {'instruction_set': instruction_set,
'assume_aligned': aligned_and_padding[0],
'nontemporal': nontemporal,
'assume_inner_stride_one': True,
'assume_sufficient_line_padding': aligned_and_padding[1]}
'assume_aligned': aligned_and_padding[0],
'nontemporal': nontemporal,
'assume_inner_stride_one': True,
'assume_sufficient_line_padding': aligned_and_padding[1]}
time_steps = 100
size1 = (64, 32)
size2 = (666, 34)
......
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