Skip to content
Snippets Groups Projects
Forked from waLBerla / waLBerla
94 commits behind the upstream repository.
04_LBComplexGeometry.py 5.04 KiB
import sympy as sp
import numpy as np
import pystencils as ps

from lbmpy.enums import SubgridScaleModel
from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Method, Stencil

from lbmpy.creationfunctions import create_lb_update_rule, create_lb_collision_rule, create_lb_method

from lbmpy.boundaries import NoSlip, FreeSlip, UBB, ExtrapolationOutflow, NoSlipLinearBouzidi, QuadraticBounceBack

from pystencils_walberla import CodeGeneration, generate_sweep, generate_info_header
from lbmpy_walberla import generate_boundary, generate_lattice_model, generate_lbm_package, lbm_boundary_generator
from lbmpy_walberla.additional_data_handler import QuadraticBounceBackAdditionalDataHandler
from pystencils.typing import TypedSymbol

with CodeGeneration() as ctx:
    data_type = "float64" if ctx.double_accuracy else "float32"
    stencil = LBStencil(Stencil.D3Q27)
    omega = sp.Symbol('omega')

    if ctx.gpu:
        target = ps.Target.GPU
    else:
        target = ps.Target.CPU

    layout = 'fzyx'

    compile_time_block_size = False

    max_threads = 256 if target == ps.Target.GPU else None

    if compile_time_block_size:
        sweep_block_size = (128, 1, 1)
    else:
        sweep_block_size = (TypedSymbol("cudaBlockSize0", np.int32),
                            TypedSymbol("cudaBlockSize1", np.int32),
                            TypedSymbol("cudaBlockSize2", np.int32))

    gpu_indexing_params = {'block_size': sweep_block_size}
    
    #   PDF Fields
    pdfs, pdfs_tmp = ps.fields(f'pdfs({stencil.Q}),pdfs_tmp({stencil.Q}): {data_type}[{stencil.D}D]', layout=layout)

    #  Output Fields
    velocity = ps.fields(f"velocity({stencil.D}):{data_type}[{stencil.D}D]", layout=layout)
    density  = ps.fields(f"density({1}):{data_type}[{stencil.D}D]", layout=layout)

    macroscopic_fields = {'density': density, 'velocity': velocity}

    # LBM Optimisation
    lbm_opt = LBMOptimisation(cse_global=True,
                              symbolic_field=pdfs,
                              symbolic_temporary_field=pdfs_tmp,
                              field_layout=layout)

    #   ==================
    #      Method Setup
    #   ==================
    lbm_config = LBMConfig(stencil=stencil,
                           method=Method.CUMULANT,
                           galilean_correction=True,
                           zero_centered=True,
                           relaxation_rate=omega,
                           fourth_order_correction=0.01,
                           compressible=True,
                           output=macroscopic_fields)


    lbm_method = create_lb_method(lbm_config=lbm_config)
    lbm_update_rule = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
    collision_rule  = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)

    freeslip = lbm_boundary_generator("FreeSlipBC", flag_uid="FreeSlip", boundary_object=FreeSlip(stencil))

    noslip = lbm_boundary_generator(class_name='NoSlipBC', flag_uid='NoSlip', boundary_object=NoSlip(),
                                    field_data_type=data_type)
    
    noslipQbb = lbm_boundary_generator(class_name='NoSlipQBBBC', flag_uid='NoSlipQBB',
                                    boundary_object=QuadraticBounceBack( relaxation_rate=omega ),
                                    field_data_type=data_type)
    
    wallDistanceCallback = lambda *args: None
    obj_noslipQbb = lbm_boundary_generator(class_name='ObjNoSlipQBBBC', flag_uid='ObjNoSlipQBB',
                                    boundary_object=QuadraticBounceBack( relaxation_rate=omega ),
                                    field_data_type=data_type)

    outflow = lbm_boundary_generator( class_name='OutflowBC', flag_uid='Outflow',
                                      boundary_object=ExtrapolationOutflow(stencil[4], lbm_method),
                                      field_data_type=data_type)

    ubb = lbm_boundary_generator(class_name='UBBBC', flag_uid='UBB',
                                 boundary_object=UBB(lambda *args: None, dim=stencil.D, data_type=data_type),
                                 field_data_type=data_type)

    generate_lbm_package(ctx, name="LBComplexGeometry",
                         collision_rule=collision_rule,
                         lbm_config=lbm_config, lbm_optimisation=lbm_opt,
                         nonuniform=False, boundaries=[outflow, ubb, freeslip, noslip, noslipQbb, obj_noslipQbb],
                         macroscopic_fields=macroscopic_fields,
                         gpu_indexing_params=gpu_indexing_params,
                         max_threads=max_threads, target=target,
                         data_type=data_type, pdfs_data_type=data_type,
                         cpu_openmp=ctx.openmp, set_pre_collision_pdfs=False)

    field_typedefs = {'VelocityField_T': velocity,
                      'ScalarField_T': density}

    stencil_typedefs = {'Stencil_T': stencil}

    # Info header containing correct template definitions for stencil and field
    generate_info_header(ctx, 'InfoHeader',stencil_typedefs=stencil_typedefs,
                         field_typedefs=field_typedefs)