Commit 58d371d7 authored by Markus Holzer's avatar Markus Holzer Committed by Christoph Schwarzmeier
Browse files

Generated outflow bc

parent c81e1923
......@@ -13,6 +13,7 @@ Dominik Bartuschat
Ehsan Fattahi
Felix Winterhalter
Florian Schornbaum
Frederik Hennig
Grigorii Drozdov
Helen Schottenhamml
Igor Ostanin
......
......@@ -69,11 +69,10 @@ options_dict = {
'relaxation_rate': omega,
},
'cumulant': {
'method': 'cumulant',
'stencil': 'D3Q19',
'compressible': True,
'method': 'mrt',
'cumulant': True,
'relaxation_rates': [0, omega, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
'relaxation_rate': omega,
},
}
......
......@@ -55,11 +55,10 @@ options_dict = {
'relaxation_rate': omega,
},
'cumulant': {
'method': 'cumulant',
'stencil': 'D3Q19',
'compressible': True,
'method': 'mrt',
'cumulant': True,
'relaxation_rates': [0, omega, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
'relaxation_rate': omega,
},
}
......
......@@ -40,11 +40,13 @@ using namespace walberla;
Field<walberla::real_t,1>,\
Field<walberla::real_t,2>,\
Field<walberla::real_t,3>,\
Field<walberla::real_t,4>,\
Field<walberla::real_t,9>,\
Field<walberla::real_t,15>,\
Field<walberla::real_t,19>,\
Field<walberla::real_t,27>,\
Field<walberla::int8_t,1>,\
Field<walberla::int32_t,1>,\
Field<walberla::int64_t,1>,\
Field<walberla::int64_t,2>,\
Field<walberla::int64_t,3>,\
......
......@@ -17,5 +17,4 @@ waLBerla_add_executable(NAME multiphaseCPU
FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp contact.cpp CalculateNormals.cpp multiphase_codegen.py
DEPENDS blockforest core field postprocessing lbm geometry timeloop gui PhaseFieldCodeGenCPU)
set_target_properties(multiphaseCPU PROPERTIES CXX_VISIBILITY_PRESET hidden)
......@@ -20,5 +20,4 @@ waLBerla_add_executable(NAME multiphaseGPU
FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp CalculateNormals.cpp contact.cu multiphase_codegen.py
DEPENDS blockforest core cuda field postprocessing lbm geometry timeloop gui PhaseFieldCodeGenGPU)
set_target_properties(multiphaseGPU PROPERTIES CXX_VISIBILITY_PRESET hidden)
......@@ -36,9 +36,8 @@ optimization = {'cse_global': True,
# ==================
lbm_params = {'stencil': stencil,
'method': 'mrt_raw',
'relaxation_rates': [0, 0, 0, omega, omega, omega, 1, 1, 1],
'cumulant': True,
'method': 'cumulant',
'relaxation_rate': omega,
'compressible': True}
lbm_update_rule = create_lb_update_rule(optimization=optimization,
......
from .boundary import generate_boundary
from .walberla_lbm_generation import RefinementScaling, generate_lattice_model
from .packinfo import generate_lb_pack_info
__all__ = ['generate_lattice_model', 'RefinementScaling', 'generate_boundary']
__all__ = ['generate_lattice_model', 'RefinementScaling', 'generate_boundary', 'generate_lb_pack_info']
from pystencils.stencil import inverse_direction
from lbmpy.advanced_streaming import AccessPdfValues, numeric_offsets, numeric_index
from lbmpy.boundaries import ExtrapolationOutflow, UBB
from pystencils_walberla.additional_data_handler import AdditionalDataHandler
class UBBAdditionalDataHandler(AdditionalDataHandler):
def __init__(self, stencil, boundary_object):
assert isinstance(boundary_object, UBB)
self._boundary_object = boundary_object
super(UBBAdditionalDataHandler, self).__init__(stencil=stencil)
@property
def constructor_arguments(self):
return ", std::function<Vector3<real_t>(const Cell &, const shared_ptr<StructuredBlockForest>&, IBlock&)>& " \
"velocityCallback "
@property
def initialiser_list(self):
return "elementInitaliser(velocityCallback),"
@property
def additional_arguments_for_fill_function(self):
return "blocks, "
@property
def additional_parameters_for_fill_function(self):
return " const shared_ptr<StructuredBlockForest> &blocks, "
def data_initialisation(self, direction):
init_list = ["Vector3<real_t> InitialisatonAdditionalData = elementInitaliser(Cell(it.x(), it.y(), it.z()), "
"blocks, *block);", "element.vel_0 = InitialisatonAdditionalData[0];",
"element.vel_1 = InitialisatonAdditionalData[1];"]
if self._dim == 3:
init_list.append("element.vel_2 = InitialisatonAdditionalData[2];")
return "\n".join(init_list)
@property
def additional_member_variable(self):
return "std::function<Vector3<real_t>(const Cell &, const shared_ptr<StructuredBlockForest>&, IBlock&)> " \
"elementInitaliser; "
class OutflowAdditionalDataHandler(AdditionalDataHandler):
def __init__(self, stencil, boundary_object, target='cpu', field_name='pdfs'):
assert isinstance(boundary_object, ExtrapolationOutflow)
self._boundary_object = boundary_object
self._stencil = boundary_object.stencil
self._lb_method = boundary_object.lb_method
self._normal_direction = boundary_object.normal_direction
self._field_name = field_name
self._target = target
super(OutflowAdditionalDataHandler, self).__init__(stencil=stencil)
assert sum([a != 0 for a in self._normal_direction]) == 1,\
"The outflow boundary is only implemented for straight walls at the moment."
@property
def constructor_arguments(self):
return f", BlockDataID {self._field_name}CPUID_" if self._target == 'gpu' else ""
@property
def initialiser_list(self):
return f"{self._field_name}CPUID({self._field_name}CPUID_)," if self._target == 'gpu' else ""
@property
def additional_field_data(self):
identifier = "CPU" if self._target == "gpu" else ""
return f"auto {self._field_name} = block->getData< field::GhostLayerField<double, " \
f"{len(self._stencil)}> >({self._field_name}{identifier}ID); "
def data_initialisation(self, direction_index):
pdf_acc = AccessPdfValues(self._boundary_object.stencil,
streaming_pattern=self._boundary_object.streaming_pattern,
timestep=self._boundary_object.zeroth_timestep,
streaming_dir='out')
init_list = []
for key, value in self.get_init_dict(pdf_acc, direction_index).items():
init_list.append(f"element.{key} = {self._field_name}->get({value});")
return "\n".join(init_list)
@property
def additional_member_variable(self):
return f"BlockDataID {self._field_name}CPUID;"
@property
def stencil_info(self):
stencil_info = []
for i, d in enumerate(self._stencil):
if any([a != 0 and b != 0 and a == b for a, b in zip(self._normal_direction, d)]):
direction = d if self._dim == 3 else d + (0,)
stencil_info.append((i, direction, ", ".join([str(e) for e in direction])))
return stencil_info
def get_init_dict(self, pdf_accessor, direction_index):
"""The Extrapolation Outflow boundary needs additional data. This function provides a list of all values
which have to be initialised"""
position = ["it.x()", "it.y()", "it.z()"]
direction = self._stencil[direction_index]
inv_dir = self._stencil.index(inverse_direction(direction))
tangential_offset = tuple(offset - normal for offset, normal in zip(direction, self._normal_direction))
result = {}
pos = []
offsets = numeric_offsets(pdf_accessor.accs[inv_dir])
for p, o, t in zip(position, offsets, tangential_offset):
pos.append(p + f" + cell_idx_c({str(o + t)})")
if self._dim == 2:
pos.append("0")
pos.append(str(numeric_index(pdf_accessor.accs[inv_dir])[0]))
result[f'pdf'] = ', '.join(pos)
pos = []
for p, o, t in zip(position, offsets, tangential_offset):
pos.append(p + f" + cell_idx_c({str(o + t)})")
if self._dim == 2:
pos.append("0")
pos.append(str(numeric_index(pdf_accessor.accs[inv_dir])[0]))
result[f'pdf_nd'] = ', '.join(pos)
return result
import pystencils_walberla.boundary
from lbmpy.boundaries.boundaryhandling import create_lattice_boltzmann_boundary_kernel
from lbmpy.advanced_streaming import Timestep, is_inplace
def generate_boundary(generation_context,
......@@ -7,15 +8,29 @@ def generate_boundary(generation_context,
boundary_object,
lb_method,
field_name='pdfs',
streaming_pattern='pull',
always_generate_separate_classes=False,
additional_data_handler=None,
**create_kernel_params):
def boundary_creation_function(field, index_field, stencil, boundary_functor, target='cpu', openmp=True, **kwargs):
return create_lattice_boltzmann_boundary_kernel(field,
index_field,
lb_method,
boundary_functor,
target=target,
**kwargs)
pargs = (field, index_field, lb_method, boundary_functor)
kwargs = {'target': target, **kwargs}
if is_inplace(streaming_pattern) or always_generate_separate_classes:
return {
'EvenSweep': create_lattice_boltzmann_boundary_kernel(*pargs,
streaming_pattern=streaming_pattern,
prev_timestep=Timestep.EVEN,
**kwargs),
'OddSweep': create_lattice_boltzmann_boundary_kernel(*pargs,
streaming_pattern=streaming_pattern,
prev_timestep=Timestep.ODD,
**kwargs)
}
else:
return create_lattice_boltzmann_boundary_kernel(*pargs,
streaming_pattern=streaming_pattern,
prev_timestep=Timestep.BOTH,
**kwargs)
pystencils_walberla.boundary.generate_boundary(generation_context,
class_name,
......@@ -25,4 +40,5 @@ def generate_boundary(generation_context,
index_shape=[len(lb_method.stencil)],
kernel_creation_function=boundary_creation_function,
namespace='lbm',
additional_data_handler=additional_data_handler,
**create_kernel_params)
from collections import defaultdict
from lbmpy.advanced_streaming.utility import Timestep, get_accessor, get_timesteps
from lbmpy.advanced_streaming.communication import _extend_dir
from pystencils.stencil import inverse_direction
from pystencils_walberla.codegen import comm_directions, generate_pack_info
from pystencils import Assignment, Field
def generate_lb_pack_info(generation_context,
class_name_prefix: str,
stencil,
pdf_field,
streaming_pattern='pull',
lb_collision_rule=None,
always_generate_separate_classes=False,
namespace='lbm',
**create_kernel_params):
"""Generates waLBerla MPI PackInfos for an LBM kernel, based on a given method
and streaming pattern. For in-place streaming patterns, two PackInfos are generated;
one for the even and another for the odd time steps.
Args:
generation_context: see documentation of `generate_sweep`
class_name_prefix: Prefix of the desired class name which will be extended with
'Even' or 'Odd' for in-place kernels
stencil: The tuple of directions specifying the employed LB stencil.
pdf_field: pdf field for which the pack info is created
streaming_pattern: The employed streaming pattern.
lb_collision_rule: Optional. The collision rule defining the lattice boltzmann kernel, as returned
by `create_lb_collision_rule`. If specified, it will be scanned for non-local
accesses to other fields other than the PDF fields (as might be required for
computing gradients in coupled simulations), whose communication will then
be included in the PackInfo.
always_generate_separate_classes: If True, generate a pair of Even/Odd PackInfos even for a two-
fields kernel (i.e. the pull/push patterns). Otherwise, for two-fields
kernels, only one PackInfo class will be generated without a
suffix to its name.
namespace: inner namespace of the generated class
**create_kernel_params: remaining keyword arguments are passed to `pystencils.create_kernel`
"""
timesteps = [Timestep.EVEN, Timestep.ODD] \
if always_generate_separate_classes \
else get_timesteps(streaming_pattern)
common_spec = defaultdict(set)
if lb_collision_rule is not None:
assignments = lb_collision_rule.all_assignments
reads = set()
for a in assignments:
if not isinstance(a, Assignment):
continue
reads.update(a.rhs.atoms(Field.Access))
for fa in reads:
assert all(abs(e) <= 1 for e in fa.offsets)
if all(offset == 0 for offset in fa.offsets):
continue
comm_direction = inverse_direction(fa.offsets)
for comm_dir in comm_directions(comm_direction):
common_spec[(comm_dir,)].add(fa.field.center(*fa.index))
for t in timesteps:
spec = common_spec.copy()
write_accesses = get_accessor(streaming_pattern, t).write(pdf_field, stencil)
for comm_dir in stencil:
if all(d == 0 for d in comm_dir):
continue
for streaming_dir in set(_extend_dir(comm_dir)) & set(stencil):
d = stencil.index(streaming_dir)
fa = write_accesses[d]
spec[(comm_dir,)].add(fa)
if t == Timestep.EVEN:
class_name_suffix = 'Even'
elif t == Timestep.ODD:
class_name_suffix = 'Odd'
else:
class_name_suffix = ''
class_name = class_name_prefix + class_name_suffix
generate_pack_info(generation_context, class_name, spec, namespace=namespace, **create_kernel_params)
......@@ -6,7 +6,6 @@ from jinja2 import Environment, PackageLoader, StrictUndefined, Template
from sympy.tensor import IndexedBase
import pystencils as ps
# from lbmpy.creationfunctions import create_lb_update_rule, update_with_default_parameters
from lbmpy.fieldaccess import CollideOnlyInplaceAccessor, StreamPullTwoFieldsAccessor
from lbmpy.relaxationrates import relaxation_rate_scaling
from lbmpy.stencils import get_stencil
......
from pystencils.stencil import inverse_direction
class AdditionalDataHandler:
"""Base class that defines how to handle boundary conditions holding additional data."""
def __init__(self, stencil):
self._dim = len(stencil[0])
# waLBerla is a 3D framework. Therefore, a zero for the z index has to be added if we work in 2D
if self._dim == 2:
self._walberla_stencil = ()
for d in stencil:
d = d + (0,)
self._walberla_stencil = self._walberla_stencil + (d,)
else:
self._walberla_stencil = stencil
@property
def constructor_arguments(self):
return ""
@property
def initialiser_list(self):
return ""
@property
def additional_arguments_for_fill_function(self):
return ""
@property
def additional_parameters_for_fill_function(self):
return ""
@property
def additional_field_data(self):
return ""
def data_initialisation(self, direction_index):
return ""
@property
def additional_member_variable(self):
return ""
@property
def stencil_info(self):
return [(i, d, ", ".join([str(e) for e in d])) for i, d in enumerate(self._walberla_stencil)]
@property
def inverse_directions(self):
inv_dirs = []
for direction in self._walberla_stencil:
inv_dirs.append(self._walberla_stencil.index(inverse_direction(direction)))
return inv_dirs
from collections import OrderedDict
import numpy as np
from jinja2 import Environment, PackageLoader, StrictUndefined
from pystencils import Field, FieldType
......@@ -8,6 +9,7 @@ from pystencils.boundaries.createindexlist import (
from pystencils.data_types import TypedSymbol, create_type
from pystencils_walberla.codegen import KernelInfo, default_create_kernel_parameters
from pystencils_walberla.jinja_filters import add_pystencils_filters_to_jinja_env
from pystencils_walberla.additional_data_handler import AdditionalDataHandler
def generate_boundary(generation_context,
......@@ -20,7 +22,12 @@ def generate_boundary(generation_context,
kernel_creation_function=None,
target='cpu',
namespace='pystencils',
additional_data_handler=None,
**create_kernel_params):
if boundary_object.additional_data and additional_data_handler is None:
raise ValueError("Boundary object has additional data but you have not provided an AdditionalDataHandler.")
struct_name = "IndexInfo"
boundary_object.name = class_name
dim = len(neighbor_stencil[0])
......@@ -44,36 +51,44 @@ def generate_boundary(generation_context,
if not kernel_creation_function:
kernel_creation_function = create_boundary_kernel
kernel = kernel_creation_function(field, index_field, neighbor_stencil, boundary_object, **create_kernel_params)
kernel.function_name = "boundary_" + boundary_object.name
kernel.assumed_inner_stride_one = False
# waLBerla is a 3D framework. Therefore, a zero for the z index has to be added if we work in 2D
if dim == 2:
stencil = ()
for d in neighbor_stencil:
d = d + (0,)
stencil = stencil + (d,)
kernels = kernel_creation_function(field, index_field, neighbor_stencil, boundary_object, **create_kernel_params)
if isinstance(kernels, dict):
sweep_to_kernel_info_dict = OrderedDict()
dummy_kernel_info = None
for sweep_class, sweep_kernel in kernels.items():
sweep_kernel.function_name = "boundary_" + boundary_object.name + '_' + sweep_class
sweep_kernel.assumed_inner_stride_one = False
kernel_info = KernelInfo(sweep_kernel)
sweep_to_kernel_info_dict[sweep_class] = kernel_info
if dummy_kernel_info is None:
dummy_kernel_info = kernel_info
# elif not dummy_kernel_info.has_same_interface(kernel_info):
# raise ValueError("Multiple boundary sweeps must have the same kernel interface!")
multi_sweep = True
else:
stencil = neighbor_stencil
multi_sweep = False
kernel = kernels
kernel.function_name = "boundary_" + boundary_object.name
kernel.assumed_inner_stride_one = False
kernel_info = KernelInfo(kernel)
sweep_to_kernel_info_dict = {'': kernel_info}
dummy_kernel_info = kernel_info
stencil_info = [(i, d, ", ".join([str(e) for e in d])) for i, d in enumerate(stencil)]
inv_dirs = []
for direction in stencil:
inverse_dir = tuple([-i for i in direction])
inv_dirs.append(stencil.index(inverse_dir))
if additional_data_handler is None:
additional_data_handler = AdditionalDataHandler(stencil=neighbor_stencil)
context = {
'class_name': boundary_object.name,
'sweep_classes': sweep_to_kernel_info_dict,
'multi_sweep': multi_sweep,
'dummy_kernel_info': dummy_kernel_info,
'StructName': struct_name,
'StructDeclaration': struct_from_numpy_dtype(struct_name, index_struct_dtype),
'kernel': KernelInfo(kernel),
'stencil_info': stencil_info,
'inverse_directions': inv_dirs,
'dim': dim,
'target': target,
'namespace': namespace,
'inner_or_boundary': boundary_object.inner_or_boundary
'inner_or_boundary': boundary_object.inner_or_boundary,
'additional_data_handler': additional_data_handler
}
env = Environment(loader=PackageLoader('pystencils_walberla'), undefined=StrictUndefined)
......@@ -83,8 +98,8 @@ def generate_boundary(generation_context,
source = env.get_template('Boundary.tmpl.cpp').render(**context)
source_extension = "cpp" if target == "cpu" else "cu"
generation_context.write_file("{}.h".format(class_name), header)
generation_context.write_file("{}.{}".format(class_name, source_extension), source)
generation_context.write_file(f"{class_name}.h", header)
generation_context.write_file(f"{class_name}.{source_extension}", source)
def generate_staggered_boundary(generation_context, class_name, boundary_object,
......@@ -102,23 +117,23 @@ def generate_staggered_flux_boundary(generation_context, class_name, boundary_ob
def struct_from_numpy_dtype(struct_name, numpy_dtype):
result = "struct %s { \n" % (struct_name,)
result = f"struct {struct_name} {{ \n"
equality_compare = []
constructor_params = []
constructor_initializer_list = []
for name, (sub_type, offset) in numpy_dtype.fields.items():
pystencils_type = create_type(sub_type)
result += " %s %s;\n" % (pystencils_type, name)
result += f" {pystencils_type} {name};\n"
if name in boundary_index_array_coordinate_names or name == direction_member_name:
constructor_params.append("%s %s_" % (pystencils_type, name))
constructor_initializer_list.append("%s(%s_)" % (name, name))
constructor_params.append(f"{pystencils_type} {name}_")
constructor_initializer_list.append(f"{name}({name}_)")
else:
constructor_initializer_list.append("%s()" % name)
constructor_initializer_list.append(f"{name}()")
if pystencils_type.is_float():
equality_compare.append("floatIsEqual(%s, o.%s)" % (name, name))
equality_compare.append(f"floatIsEqual({name}, o.{name})")
else:
equality_compare.append("%s == o.%s" % (name, name))
equality_compare.append(f"{name} == o.{name}")
result += " %s(%s) : %s {}\n" % \
(struct_name, ", ".join(constructor_params), ", ".join(constructor_initializer_list))
......
......@@ -20,7 +20,7 @@ __all__ = ['generate_sweep', 'generate_pack_info', 'generate_pack_info_for_field
def generate_sweep(generation_context, class_name, assignments,
namespace='pystencils', field_swaps=(), staggered=False, varying_parameters=(),
inner_outer_split=False,
inner_outer_split=False, ghost_layers_to_include=0,
**create_kernel_params):
"""Generates a waLBerla sweep from a pystencils representation.
......@@ -44,6 +44,8 @@ def generate_sweep(generation_context, class_name, assignments,
the C++ class constructor even if the kernel does not need them.
inner_outer_split: if True generate a sweep that supports separate iteration over inner and outer regions
to allow for communication hiding.
ghost_layers_to_include: determines how many ghost layers should be included for the Sweep.
This is relevant if a setter kernel should also set correct values to the ghost layers.
**create_kernel_params: remaining keyword arguments are passed to `pystencils.create_kernel`
"""
create_kernel_params = default_create_kernel_parameters(generation_context, create_kernel_params)
......@@ -78,6 +80,7 @@ def generate_sweep(generation_context, class_name, assignments,
'class_name': class_name,
'target': create_kernel_params.get("target", "cpu"),
'headers': get_headers(ast),
'ghost_layers_to_include': ghost_layers_to_include
}
header = env.get_template("Sweep.tmpl.h").render(**jinja_context)
source = env.get_template("Sweep.tmpl.cpp").render(**jinja_context)
......@@ -93,6 +96,7 @@ def generate_sweep(generation_context, class_name, assignments,
'target': create_kernel_params.get("target", "cpu"),
'field': representative_field,
'headers': get_headers(ast),
'ghost_layers_to_include': 0
}
header = env.get_template("SweepInnerOuter.tmpl.h").render(**jinja_context)
source = env.get_template("SweepInnerOuter.tmpl.cpp").render(**jinja_context)
......@@ -132,7 +136,7 @@ def generate_pack_info_from_kernel(generation_context, class_name: str, assignme
class_name: name of the generated class
assignments: list of assignments from the compute kernel - generates PackInfo for "pull" part only
i.e. the kernel is expected to only write to the center
kind:
kind: can either be pull or push
**create_kernel_params: remaining keyword arguments are passed to `pystencils.create_kernel`
"""
assert kind in ('push', 'pull')
......@@ -240,7 +244,6 @@ def generate_pack_info(generation_context, class_name: str,
pack_kernels[direction_strings] = KernelInfo(pack_ast)
unpack_kernels[direction_strings] = KernelInfo(unpack_ast)