Commit f59f899e authored by Sebastian Bindgen 's avatar Sebastian Bindgen Committed by Markus Holzer
Browse files

Revert "Use regular periodic boundaries"

This reverts commit 38d9aa5a.
parent 18670f8f
......@@ -13,6 +13,8 @@ tests-and-coverage:
- $ENABLE_NIGHTLY_BUILDS
image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
script:
- env
- pip list
- export NUM_CORES=$(nproc --all)
- mkdir -p ~/.config/matplotlib
- echo "backend:template" > ~/.config/matplotlib/matplotlibrc
......@@ -38,6 +40,8 @@ test-longrun:
- $ENABLE_NIGHTLY_BUILDS
image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
script:
- env
- pip list
- export NUM_CORES=$(nproc --all)
- mkdir -p ~/.config/matplotlib
- echo "backend:template" > ~/.config/matplotlib/matplotlibrc
......@@ -103,6 +107,8 @@ pycodegen-integration:
when: manual
allow_failure: true
script:
- env
- pip list
- git clone https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pycodegen.git
- cd pycodegen
- git submodule sync --recursive
......
This diff is collapsed.
......@@ -21,4 +21,4 @@ You can open the notebooks directly to play around with the code examples.
/notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb
/notebooks/demo_automatic_chapman_enskog_analysis.ipynb
/notebooks/demo_thermalized_lbm.ipynb
/notebooks/demo_theoretical_background_generic_equilibrium_construction.ipynb
\ No newline at end of file
/notebooks/demo_theoretical_background_generic_equilibrium_construction.ipynb
This source diff could not be displayed because it is too large. You can view the blob instead.
from functools import partial
import pystencils as ps
from pystencils.astnodes import LoopOverCoordinate
from pystencils.slicing import get_periodic_boundary_functor
from lbmpy.creationfunctions import create_lb_update_rule
from lbmpy.stencils import get_stencil
from lbmpy.updatekernels import create_stream_pull_with_output_kernel
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
from lbmpy.relaxationrates import lattice_viscosity_from_relaxation_rate
import sympy as sp
import numpy as np
def get_le_boundary_functor(neighbor_stencil, shear_offset, ghost_layers=1, thickness=None, n=64):
functor_2 = get_periodic_boundary_functor(neighbor_stencil, ghost_layers, thickness)
def functor(pdfs, **_):
functor_2(pdfs)
weight = np.fmod(shear_offset[0] + n, 1.)
# First, we interpolate the upper pdfs
for i in range(len(pdfs[:, ghost_layers, :])):
ind1 = int(np.floor(i - shear_offset[0]) % n)
ind2 = int(np.ceil(i - shear_offset[0]) % n)
res = (1 - weight) * pdfs[ind1, ghost_layers, :] + weight * pdfs[ind2, ghost_layers, :]
pdfs[i, -ghost_layers, :] = res
# Second, we interpolate the lower pdfs
for i in range(len(pdfs[:, -ghost_layers, :])):
ind1 = int(np.floor(i + shear_offset[0]) % n)
ind2 = int(np.ceil(i + shear_offset[0]) % n)
res = (1 - weight) * pdfs[ind1, -ghost_layers - 1, :] + weight * pdfs[ind2, -ghost_layers - 1, :]
pdfs[i, ghost_layers - 1, :] = res
return functor
def get_solution_navier_stokes(x, t, viscosity, velocity=1.0, height=1.0, max_iterations=10):
w = x / height - 0.5
for k in np.arange(1, max_iterations + 1):
w += 1.0 / (np.pi * k) * np.exp(-4 * np.pi ** 2 * viscosity * k ** 2 / height ** 2 * t) * np.sin(
2 * np.pi / height * k * x)
return velocity * w
def test_lees_edwards():
domain_size = (64, 64)
omega = 1.0 # relaxation rate of first component
shear_velocity = 0.1 # shear velocity
shear_dir = 0 # direction of shear flow
shear_dir_normal = 1 # direction normal to shear plane, for interpolation
stencil = get_stencil("D2Q9")
dim = len(stencil[0])
dh = ps.create_data_handling(domain_size, periodicity=True, default_target='cpu')
src = dh.add_array('src', values_per_cell=len(stencil))
dh.fill('src', 1.0, ghost_layers=True)
dst = dh.add_array_like('dst', 'src')
dh.fill('dst', 0.0, ghost_layers=True)
force = dh.add_array('force', values_per_cell=dh.dim)
dh.fill('force', 0.0, ghost_layers=True)
rho = dh.add_array('rho', values_per_cell=1)
dh.fill('rho', 1.0, ghost_layers=True)
u = dh.add_array('u', values_per_cell=dh.dim)
dh.fill('u', 0.0, ghost_layers=True)
counters = [LoopOverCoordinate.get_loop_counter_symbol(i) for i in range(dim)]
points_up = sp.Symbol('points_up')
points_down = sp.Symbol('points_down')
u_p = sp.Piecewise((1, sp.And(ps.data_types.type_all_numbers(counters[1] <= 1, 'int'), points_down)),
(-1, sp.And(ps.data_types.type_all_numbers(counters[1] >= src.shape[1] - 2, 'int'),
points_up)), (0, True)) * shear_velocity
collision = create_lb_update_rule(stencil=stencil,
relaxation_rate=omega,
compressible=True,
velocity_input=u.center_vector + sp.Matrix([u_p, 0]),
density_input=rho,
force_model='luo',
force=force.center_vector,
kernel_type='collide_only',
optimization={'symbolic_field': src})
to_insert = [s.lhs for s in collision.subexpressions
if collision.method.first_order_equilibrium_moment_symbols[shear_dir]
in s.free_symbols]
for s in to_insert:
collision = collision.new_with_inserted_subexpression(s)
ma = []
for a, c in zip(collision.main_assignments, collision.method.stencil):
if c[shear_dir_normal] == -1:
b = (True, False)
elif c[shear_dir_normal] == 1:
b = (False, True)
else:
b = (False, False)
a = ps.Assignment(a.lhs, a.rhs.replace(points_down, b[0]))
a = ps.Assignment(a.lhs, a.rhs.replace(points_up, b[1]))
ma.append(a)
collision.main_assignments = ma
stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
{'density': rho, 'velocity': u})
stream_kernel = ps.create_kernel(stream, target=dh.default_target).compile()
collision_kernel = ps.create_kernel(collision, target=dh.default_target).compile()
init = macroscopic_values_setter(collision.method, velocity=(0, 0),
pdfs=src.center_vector, density=rho.center)
init_kernel = ps.create_kernel(init, ghost_layers=0).compile()
offset = [0.0]
sync_pdfs = dh.synchronization_function([src.name],
functor=partial(get_le_boundary_functor, shear_offset=offset))
dh.run_kernel(init_kernel)
time = 500
dh.all_to_gpu()
for i in range(time):
dh.run_kernel(collision_kernel)
sync_pdfs()
dh.run_kernel(stream_kernel)
dh.swap(src.name, dst.name)
offset[0] += shear_velocity
dh.all_to_cpu()
nu = lattice_viscosity_from_relaxation_rate(omega)
h = domain_size[0]
k_max = 100
analytical_solution = get_solution_navier_stokes(np.linspace(0.5, h - 0.5, h), time, nu, shear_velocity, h, k_max)
np.testing.assert_array_almost_equal(analytical_solution, dh.gather_array(u.name)[0, :, 0], decimal=5)
dh.fill(rho.name, 1.0, ghost_layers=True)
dh.run_kernel(init_kernel)
dh.fill(u.name, 0.0, ghost_layers=True)
dh.fill('force', 0.0, ghost_layers=True)
dh.cpu_arrays[force.name][64 // 3, 1, :] = [1e-2, -1e-1]
offset[0] = 0
time = 20
dh.all_to_gpu()
for i in range(time):
dh.run_kernel(collision_kernel)
sync_pdfs()
dh.run_kernel(stream_kernel)
dh.swap(src.name, dst.name)
dh.all_to_cpu()
vel_unshifted = np.array(dh.gather_array(u.name)[:, -3:-1, :])
dh.fill(rho.name, 1.0, ghost_layers=True)
dh.run_kernel(init_kernel)
dh.fill(u.name, 0.0, ghost_layers=True)
dh.fill('force', 0.0, ghost_layers=True)
dh.cpu_arrays[force.name][64 // 3, 1, :] = [1e-2, -1e-1]
offset[0] = 10
time = 20
dh.all_to_gpu()
for i in range(time):
dh.run_kernel(collision_kernel)
sync_pdfs()
dh.run_kernel(stream_kernel)
dh.swap(src.name, dst.name)
dh.all_to_cpu()
vel_shifted = np.array(dh.gather_array(u.name)[:, -3:-1, :])
vel_rolled = np.roll(vel_shifted, -offset[0], axis=0)
np.testing.assert_array_almost_equal(vel_unshifted, vel_rolled)
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -793,7 +793,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
"version": "3.8.2"
}
},
"nbformat": 4,
......
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