Commit 6e3baa97 authored by Christoph Schwarzmeier's avatar Christoph Schwarzmeier
Browse files

Merge branch 'codegen-complex-geometry-freeslip' into 'master'

FreeSlip in Codegen for complex geometry

See merge request walberla/walberla!520
parents 6266cf80 7cc9f764
Pipeline #38147 passed with stages
in 156 minutes and 54 seconds
from pystencils import Target
from pystencils.stencil import inverse_direction
from pystencils.stencil import inverse_direction, offset_to_direction_string
from lbmpy.advanced_streaming import AccessPdfValues, numeric_offsets, numeric_index
from lbmpy.advanced_streaming.indexing import MirroredStencilDirections
from lbmpy.boundaries.boundaryconditions import LbBoundary
from lbmpy.boundaries import ExtrapolationOutflow, UBB
from lbmpy.boundaries import ExtrapolationOutflow, FreeSlip, UBB
from pystencils_walberla.additional_data_handler import AdditionalDataHandler
......@@ -12,7 +14,9 @@ def default_additional_data_handler(boundary_obj: LbBoundary, lb_method, field_n
if not boundary_obj.additional_data:
return None
if isinstance(boundary_obj, UBB):
if isinstance(boundary_obj, FreeSlip):
return FreeSlipAdditionalDataHandler(lb_method.stencil, boundary_obj)
elif isinstance(boundary_obj, UBB):
return UBBAdditionalDataHandler(lb_method.stencil, boundary_obj)
elif isinstance(boundary_obj, ExtrapolationOutflow):
return OutflowAdditionalDataHandler(lb_method.stencil, boundary_obj, target=target, field_name=field_name)
......@@ -20,6 +24,90 @@ def default_additional_data_handler(boundary_obj: LbBoundary, lb_method, field_n
raise ValueError(f"No default AdditionalDataHandler available for boundary of type {boundary_obj.__class__}")
class FreeSlipAdditionalDataHandler(AdditionalDataHandler):
def __init__(self, stencil, boundary_object):
assert isinstance(boundary_object, FreeSlip)
self._boundary_object = boundary_object
super(FreeSlipAdditionalDataHandler, self).__init__(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 ""
def data_initialisation(self, direction):
def array_pattern(dtype, name, content):
return f"const {str(dtype)} {name} [] = {{ {','.join(str(c) for c in content)} }};"
offset = self._walberla_stencil[direction]
inv_offset = inverse_direction(self._walberla_stencil[direction])
offset_dtype = "int32_t"
mirror_stencil = MirroredStencilDirections.mirror_stencil
axis_mirrored = []
for i, name in enumerate(["x", "y", "z"]):
mirrored_dir = [self._walberla_stencil.index(mirror_stencil(d, i)) for d in self._walberla_stencil]
axis_mirrored.append(array_pattern(offset_dtype, f"{name}_axis_mirrored_stencil_dir", mirrored_dir))
init_list = axis_mirrored[0:self._dim]
init_list += [f"const Cell n = it.cell() + Cell({offset[0]}, {offset[1]}, {offset[2]});",
f"int32_t ref_dir = {self._walberla_stencil.index(inv_offset)}; // dir: {direction}",
"element.wnx = 0; // compute discrete normal vector of free slip wall",
"element.wny = 0;",
f"if( flagField->isPartOfMaskSet( n.x() + {inv_offset[0]}, n.y(), n.z(), domainFlag ) )",
"{",
f" element.wnx = {inv_offset[0]};",
" ref_dir = x_axis_mirrored_stencil_dir[ ref_dir ];",
"}",
f"if( flagField->isPartOfMaskSet( n.x(), n.y() + {inv_offset[1]}, n.z(), domainFlag ) )",
"{",
f" element.wny = {inv_offset[1]};",
" ref_dir = y_axis_mirrored_stencil_dir[ ref_dir ];",
"}"]
if self._dim == 3:
init_list += ["element.wnz = 0;",
f"if( flagField->isPartOfMaskSet( n.x(), n.y(), n.z() + {inv_offset[2]}, domainFlag ) )",
"{",
f" element.wnz = {inv_offset[2]};",
" ref_dir = z_axis_mirrored_stencil_dir[ ref_dir ];",
"}",
"// concave corner (neighbors are non-fluid)",
"if( element.wnx == 0 && element.wny == 0 && element.wnz == 0 )",
"{",
f" element.wnx = {inv_offset[0]};",
f" element.wny = {inv_offset[1]};",
f" element.wnz = {inv_offset[2]};",
f" ref_dir = {direction};",
"}"]
elif self._dim == 2:
init_list += ["// concave corner (neighbors are non-fluid)",
"if( element.wnx == 0 && element.wny == 0 )",
"{",
f" element.wnx = {inv_offset[0]};",
f" element.wny = {inv_offset[1]};",
f" ref_dir = {direction};",
"}"]
init_list.append("element.ref_dir = ref_dir;")
return "\n".join(init_list)
@property
def additional_member_variable(self):
return ""
class UBBAdditionalDataHandler(AdditionalDataHandler):
def __init__(self, stencil, boundary_object):
assert isinstance(boundary_object, UBB)
......
......@@ -91,6 +91,18 @@ waLBerla_generate_target_from_python(NAME GeneratedOutflowBCGenerated
waLBerla_compile_test( FILES codegen/GeneratedOutflowBC.cpp DEPENDS GeneratedOutflowBCGenerated)
waLBerla_execute_test( NAME GeneratedOutflowBC COMMAND $<TARGET_FILE:GeneratedOutflowBC> ${CMAKE_CURRENT_SOURCE_DIR}/codegen/GeneratedOutflowBC.prm )
waLBerla_generate_target_from_python(NAME GeneratedFreeSlipGenerated
FILE codegen/GeneratedFreeSlip.py
OUT_FILES GeneratedFreeSlip_Sweep.cpp GeneratedFreeSlip_Sweep.h
GeneratedFreeSlip_MacroSetter.cpp GeneratedFreeSlip_MacroSetter.h
GeneratedFreeSlip_NoSlip.cpp GeneratedFreeSlip_NoSlip.h
GeneratedFreeSlip_FreeSlip.cpp GeneratedFreeSlip_FreeSlip.h
GeneratedFreeSlip_PackInfoEven.cpp GeneratedFreeSlip_PackInfoEven.h
GeneratedFreeSlip_PackInfoOdd.cpp GeneratedFreeSlip_PackInfoOdd.h
GeneratedFreeSlip.h)
waLBerla_compile_test( FILES codegen/GeneratedFreeSlip.cpp DEPENDS GeneratedFreeSlipGenerated)
waLBerla_execute_test( NAME GeneratedFreeSlip COMMAND $<TARGET_FILE:GeneratedFreeSlip> ${CMAKE_CURRENT_SOURCE_DIR}/codegen/GeneratedFreeSlip.prm )
waLBerla_generate_target_from_python(NAME LbCodeGenerationExampleGenerated
FILE codegen/LbCodeGenerationExample.py
......
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla is free software: you can
// redistribute it and/or modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation, either version 3 of
// the License, or (at your option) any later version.
//
// waLBerla is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License along
// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file GeneratedFreeSlip.cpp
//! \author Markus Holzer <markus.holzer@fau.de>
//! \brief Periodic flow in x direction with FreeSlip in Pipe setup. Checks if there are no gradients in the
//! Velocity field in the end. Works for all streaming patterns (pull, push, aa, esotwist)
//
//======================================================================================================================
#include "blockforest/Initialization.h"
#include "blockforest/communication/UniformBufferedScheme.h"
#include "core/Environment.h"
#include "core/timing/RemainingTimeLogger.h"
#include "field/AddToStorage.h"
#include "field/vtk/VTKWriter.h"
#include "geometry/InitBoundaryHandling.h"
#include "timeloop/SweepTimeloop.h"
// Generated Files
#include "GeneratedFreeSlip.h"
using namespace walberla;
using PackInfoEven_T = lbm::GeneratedFreeSlip_PackInfoEven;
using PackInfoOdd_T = lbm::GeneratedFreeSlip_PackInfoOdd;
using flag_t = walberla::uint8_t;
using FlagField_T = FlagField< flag_t >;
auto pdfFieldAdder = [](IBlock* const block, StructuredBlockStorage* const storage) {
return new PdfField_T(storage->getNumberOfXCells(*block), storage->getNumberOfYCells(*block),
storage->getNumberOfZCells(*block), uint_t(1), field::fzyx,
make_shared< field::AllocateAligned< real_t, 64 > >());
};
//////////
// MAIN //
//////////
class AlternatingBeforeFunction
{
public:
typedef std::function< void() > BeforeFunction;
AlternatingBeforeFunction(BeforeFunction evenFunc, BeforeFunction oddFunc,
std::shared_ptr< lbm::TimestepTracker >& tracker)
: tracker_(tracker), funcs_{ evenFunc, oddFunc } {};
void operator()() { funcs_[tracker_->getCounter()](); }
private:
std::shared_ptr< lbm::TimestepTracker > tracker_;
std::vector< BeforeFunction > funcs_;
};
int main(int argc, char** argv)
{
walberla::Environment walberlaEnv(argc, argv);
auto blocks = blockforest::createUniformBlockGridFromConfig(walberlaEnv.config());
auto domainSetup = walberlaEnv.config()->getOneBlock("DomainSetup");
Vector3< uint_t > cellsPerBlock = domainSetup.getParameter< Vector3< uint_t > >("cellsPerBlock");
// read parameters
auto parameters = walberlaEnv.config()->getOneBlock("Parameters");
const real_t omega = parameters.getParameter< real_t >("omega", real_c(1.4));
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
const double remainingTimeLoggerFrequency =
parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
// create fields
BlockDataID pdfFieldID = blocks->addStructuredBlockData< PdfField_T >(pdfFieldAdder, "PDFs");
BlockDataID velFieldID = field::addToStorage< VelocityField_T >(blocks, "velocity", real_t(0), field::fzyx);
BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(1.0), field::fzyx);
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
pystencils::GeneratedFreeSlip_MacroSetter setterSweep(densityFieldID, pdfFieldID, velFieldID);
for (auto& block : *blocks)
setterSweep(&block);
// create and initialize boundary handling
auto SpecialSetups = walberlaEnv.config()->getOneBlock("SpecialSetups");
bool tubeSetup = SpecialSetups.getParameter< bool >("tubeSetup", false);
bool slopedWall = SpecialSetups.getParameter< bool >("slopedWall", false);
const FlagUID fluidFlagUID("Fluid");
const FlagUID wallFlagUID("FreeSlip");
auto boundariesConfig = walberlaEnv.config()->getOneBlock("Boundaries");
lbm::GeneratedFreeSlip_NoSlip noSlip(blocks, pdfFieldID);
lbm::GeneratedFreeSlip_FreeSlip freeSlip(blocks, pdfFieldID);
geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldId, boundariesConfig);
if (tubeSetup || slopedWall)
{
for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
FlagField_T* flagField = blockIt->template getData< FlagField_T >(flagFieldId);
auto wallFlag = flagField->getOrRegisterFlag(wallFlagUID);
auto My = blocks->getDomainCellBB().yMax() / 2.0;
auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, {
Cell globalCell = Cell(x, y, z);
blocks->transformBlockLocalToGlobalCell(globalCell, *blockIt);
if (tubeSetup)
{
real_t R = sqrt((globalCell[1] - My) * (globalCell[1] - My) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
if (R > real_c(cellsPerBlock[1]) * real_c(0.5)) addFlag(flagField->get(x, y, z), wallFlag);
}
if (slopedWall)
{
if (globalCell[2] - globalCell[1] > 0)
{
addFlag(flagField->get(x, y, z), wallFlag);
}
}
}) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
}
}
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID);
noSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("NoSlip"), fluidFlagUID);
freeSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("FreeSlip"), fluidFlagUID);
// create time loop
SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
// create communication for PdfField
blockforest::communication::UniformBufferedScheme< Stencil_T > evenComm(blocks);
evenComm.addPackInfo(make_shared< PackInfoEven_T >(pdfFieldID));
blockforest::communication::UniformBufferedScheme< Stencil_T > oddComm(blocks);
oddComm.addPackInfo(make_shared< PackInfoOdd_T >(pdfFieldID));
// Timestep Tracking and Sweeps
auto tracker = make_shared< lbm::TimestepTracker >(0);
AlternatingBeforeFunction communication(evenComm, oddComm, tracker);
lbm::GeneratedFreeSlip_Sweep UpdateSweep(densityFieldID, pdfFieldID, velFieldID, omega);
// add LBM sweep and communication to time loop
timeloop.add() << BeforeFunction(communication, "communication") << Sweep(noSlip.getSweep(tracker), "noSlip boundary");
timeloop.add() << Sweep(freeSlip.getSweep(tracker), "freeSlip boundary");
timeloop.add() << BeforeFunction(tracker->getAdvancementFunction(), "Timestep Advancement")
<< Sweep(UpdateSweep.getSweep(tracker), "LB update rule");
// log remaining time
timeloop.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
"remaining time logger");
// VTK Writer
uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0);
if (vtkWriteFrequency > 0)
{
auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "GeneratedFreeSlip_VTK", vtkWriteFrequency, 1, false,
"vtk_out", "simulation_step", false, true, true, false, 0);
auto velWriter = make_shared< field::VTKWriter< VelocityField_T > >(velFieldID, "velocity");
auto densityWriter = make_shared< field::VTKWriter< ScalarField_T > >(densityFieldID, "density");
auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(flagFieldId, "flag");
vtkOutput->addCellDataWriter(flagWriter);
vtkOutput->addCellDataWriter(velWriter);
vtkOutput->addCellDataWriter(densityWriter);
timeloop.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
}
timeloop.run();
// simple check is only valid for tube
if (tubeSetup)
{
for (auto& block : *blocks)
{
auto velField = block.getData< VelocityField_T >(velFieldID);
for (cell_idx_t i = 0; i < cell_idx_c(cellsPerBlock[1] - 1); ++i)
{
real_t v1 = velField->get(cell_idx_c(cellsPerBlock[0] / 2), i, cell_idx_c(cellsPerBlock[2] / 2), 0);
real_t v2 = velField->get(cell_idx_c(cellsPerBlock[0] / 2), i + 1, cell_idx_c(cellsPerBlock[2] / 2), 0);
real_t grad = v2 - v1;
// WALBERLA_LOG_DEVEL_VAR(grad)
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(grad, 0.0, 1e-16)
}
}
}
return EXIT_SUCCESS;
}
Parameters
{
omega 1.8;
timesteps 2001;
tubeSetup true;
vtkWriteFrequency 100;
remainingTimeLoggerFrequency 3; // in seconds
}
SpecialSetups
{
tubeSetup true;
slopedWall false;
}
DomainSetup
{
blocks < 1, 1, 1 >;
cellsPerBlock < 20, 10, 10 >;
periodic < 1, 0, 0 >;
}
Boundaries
{
Border { direction N; walldistance -1; flag FreeSlip; }
Border { direction S; walldistance -1; flag FreeSlip; }
Border { direction B; walldistance -1; flag FreeSlip; }
Border { direction T; walldistance -1; flag FreeSlip; }
}
from dataclasses import replace
from pystencils.field import fields
from lbmpy.advanced_streaming.utility import get_timesteps, is_inplace
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Stencil, Method
from lbmpy.creationfunctions import create_lb_method, create_lb_collision_rule
from lbmpy.boundaries import NoSlip, FreeSlip
from lbmpy_walberla.additional_data_handler import FreeSlipAdditionalDataHandler
from pystencils_walberla import CodeGeneration, generate_sweep, generate_info_header
from lbmpy_walberla import generate_boundary, generate_lb_pack_info, generate_alternating_lbm_boundary, generate_alternating_lbm_sweep
import sympy as sp
stencil = LBStencil(Stencil.D3Q19)
pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): double[{stencil.D}D]", layout='fzyx')
velocity_field, density_field = fields(f"velocity({stencil.D}), density(1) : double[{stencil.D}D]", layout='fzyx')
omega = sp.Symbol("omega")
output = {
'density': density_field,
'velocity': velocity_field
}
streaming_pattern = 'esotwist'
timesteps = get_timesteps(streaming_pattern)
lbm_config = LBMConfig(method=Method.SRT, stencil=stencil, relaxation_rate=omega, force=(1e-5, 0, 0),
output=output, streaming_pattern=streaming_pattern)
lbm_opt = LBMOptimisation(symbolic_field=pdfs,
cse_global=False, cse_pdfs=False)
method = create_lb_method(lbm_config=lbm_config)
# getter & setter
setter_assignments = macroscopic_values_setter(method, density=density_field.center,
velocity=velocity_field.center_vector,
pdfs=pdfs, streaming_pattern=streaming_pattern,
previous_timestep=timesteps[0])
if not is_inplace(streaming_pattern):
lbm_opt = replace(lbm_opt, symbolic_temporary_field=pdfs_tmp)
field_swaps = [(pdfs, pdfs_tmp)]
else:
field_swaps = []
collision_rule = create_lb_collision_rule(lb_method=method, lbm_config=lbm_config, lbm_optimisation=lbm_opt)
stencil_typedefs = {'Stencil_T': stencil}
field_typedefs = {'PdfField_T': pdfs,
'VelocityField_T': velocity_field,
'ScalarField_T': density_field}
with CodeGeneration() as ctx:
# sweeps
generate_alternating_lbm_sweep(ctx, 'GeneratedFreeSlip_Sweep', collision_rule, lbm_config=lbm_config,
lbm_optimisation=lbm_opt, field_swaps=field_swaps)
generate_sweep(ctx, 'GeneratedFreeSlip_MacroSetter', setter_assignments)
# boundaries
generate_alternating_lbm_boundary(ctx, 'GeneratedFreeSlip_NoSlip', NoSlip(), method,
streaming_pattern=streaming_pattern)
free_slip = FreeSlip(stencil=stencil)
free_slip_data_handler = FreeSlipAdditionalDataHandler(stencil, free_slip)
generate_alternating_lbm_boundary(ctx, 'GeneratedFreeSlip_FreeSlip', free_slip, method,
additional_data_handler=free_slip_data_handler,
streaming_pattern=streaming_pattern)
# communication
generate_lb_pack_info(ctx, 'GeneratedFreeSlip_PackInfo', stencil, pdfs, streaming_pattern=streaming_pattern,
always_generate_separate_classes=True)
# Info header containing correct template definitions for stencil and field
generate_info_header(ctx, "GeneratedFreeSlip.h",
stencil_typedefs=stencil_typedefs, field_typedefs=field_typedefs)
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