### Revert "Use regular periodic boundaries"

`This reverts commit 38d9aa5a4ace8271a2ae3c94d4ab5c4fea0120bb.`
parent 18670f8f
 ... @@ -13,6 +13,8 @@ tests-and-coverage: ... @@ -13,6 +13,8 @@ tests-and-coverage: - \$ENABLE_NIGHTLY_BUILDS - \$ENABLE_NIGHTLY_BUILDS image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full script: script: - env - pip list - export NUM_CORES=\$(nproc --all) - export NUM_CORES=\$(nproc --all) - mkdir -p ~/.config/matplotlib - mkdir -p ~/.config/matplotlib - echo "backend:template" > ~/.config/matplotlib/matplotlibrc - echo "backend:template" > ~/.config/matplotlib/matplotlibrc ... @@ -38,6 +40,8 @@ test-longrun: ... @@ -38,6 +40,8 @@ test-longrun: - \$ENABLE_NIGHTLY_BUILDS - \$ENABLE_NIGHTLY_BUILDS image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full script: script: - env - pip list - export NUM_CORES=\$(nproc --all) - export NUM_CORES=\$(nproc --all) - mkdir -p ~/.config/matplotlib - mkdir -p ~/.config/matplotlib - echo "backend:template" > ~/.config/matplotlib/matplotlibrc - echo "backend:template" > ~/.config/matplotlib/matplotlibrc ... @@ -103,6 +107,8 @@ pycodegen-integration: ... @@ -103,6 +107,8 @@ pycodegen-integration: when: manual when: manual allow_failure: true allow_failure: true script: script: - env - pip list - git clone https://gitlab-ci-token:\${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pycodegen.git - git clone https://gitlab-ci-token:\${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pycodegen.git - cd pycodegen - cd pycodegen - git submodule sync --recursive - git submodule sync --recursive ... ...
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
 ... @@ -21,4 +21,4 @@ You can open the notebooks directly to play around with the code examples. ... @@ -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_moments_cumulants_and_maxwellian_equilibrium.ipynb /notebooks/demo_automatic_chapman_enskog_analysis.ipynb /notebooks/demo_automatic_chapman_enskog_analysis.ipynb /notebooks/demo_thermalized_lbm.ipynb /notebooks/demo_thermalized_lbm.ipynb /notebooks/demo_theoretical_background_generic_equilibrium_construction.ipynb /notebooks/demo_theoretical_background_generic_equilibrium_construction.ipynb \ No newline at end of file
This diff is collapsed.
This diff is collapsed.
 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 + n, 1.) # First, we interpolate the upper pdfs for i in range(len(pdfs[:, ghost_layers, :])): ind1 = int(np.floor(i - shear_offset) % n) ind2 = int(np.ceil(i - shear_offset) % 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) % n) ind2 = int(np.ceil(i + shear_offset) % 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) 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, 'int'), points_down)), (-1, sp.And(ps.data_types.type_all_numbers(counters >= src.shape - 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)) a = ps.Assignment(a.lhs, a.rhs.replace(points_up, b)) 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 += shear_velocity dh.all_to_cpu() nu = lattice_viscosity_from_relaxation_rate(omega) h = domain_size 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 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 = 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, axis=0) np.testing.assert_array_almost_equal(vel_unshifted, vel_rolled)
This diff is collapsed.
 ... @@ -793,7 +793,7 @@ ... @@ -793,7 +793,7 @@ "name": "python", "name": "python", "nbconvert_exporter": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "pygments_lexer": "ipython3", "version": "3.6.8" "version": "3.8.2" } } }, }, "nbformat": 4, "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!
{"ft_lang_label":"__label__en","ft_lang_prob":0.92750025,"math_prob":0.8673384,"size":378,"snap":"2022-05-2022-21","text_gpt3_token_len":105,"char_repetition_ratio":0.24598931,"word_repetition_ratio":0.40816328,"special_character_ratio":0.24338624,"punctuation_ratio":0.11666667,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.98394865,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2022-05-19T02:43:39Z\",\"WARC-Record-ID\":\"<urn:uuid:7dfa29b4-c53f-4a10-a8b6-65bd7e8c2c72>\",\"Content-Length\":\"544238\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:8bfd2d78-b0bb-4896-ac5c-6e9c993c399c>\",\"WARC-Concurrent-To\":\"<urn:uuid:726e91af-5bac-447f-9772-3e5766717d7f>\",\"WARC-IP-Address\":\"131.188.39.6\",\"WARC-Target-URI\":\"https://i10git.cs.fau.de/pycodegen/lbmpy/-/commit/f59f899eb34406a695f41c3b6fd39d3cee394ca0?view=parallel\",\"WARC-Payload-Digest\":\"sha1:I7L3HPUZVAZ7PQ7WBEHRECQ2MZGQI5X2\",\"WARC-Block-Digest\":\"sha1:GIV4WO2IDOJAMF4654WKRMG6KOM4SGP7\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2022/CC-MAIN-2022-21/CC-MAIN-2022-21_segments_1652662522741.25_warc_CC-MAIN-20220519010618-20220519040618-00000.warc.gz\"}"}