Skip to content
Snippets Groups Projects
Commit d514fba4 authored by Markus Holzer's avatar Markus Holzer Committed by Philipp Suffa
Browse files

Integration for Interpolation BCs

parent 50a01d70
Branches
Tags
No related merge requests found
......@@ -9,21 +9,38 @@ try:
except ImportError:
from lbmpy.custom_code_nodes import MirroredStencilDirections
from lbmpy.boundaries.boundaryconditions import LbBoundary
from lbmpy.boundaries import ExtrapolationOutflow, FreeSlip, UBB, DiffusionDirichlet
from lbmpy.boundaries import (ExtrapolationOutflow, FreeSlip, UBB, DiffusionDirichlet,
NoSlipLinearBouzidi, QuadraticBounceBack)
from pystencils_walberla.additional_data_handler import AdditionalDataHandler
def default_additional_data_handler(boundary_obj: LbBoundary, lb_method, field_name, target=Target.CPU):
interpolation_bc_check_template = """
if(!isFlagSet(it.neighbor({cx}, {cy}, {cz}, 0), domainFlag)){{
//Linear-Bouzidi requires 2 fluid nodes: if the 2nd node is not available abort,
//apply Bounce Back at that point. This clearly lowers the accuracy and makes inconsistent the
//calculation of the total force
element.q = -1.0;
WALBERLA_LOG_INFO_ON_ROOT("Warning: Bouzidi cannot be applied at least on one boundary link.")
}} //end if to check Bouzidi applicability
"""
def default_additional_data_handler(boundary_obj: LbBoundary, lb_method, field_name, target=Target.CPU,
pdfs_data_type=None, zeroth_timestep=None):
if not boundary_obj.additional_data:
return None
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)
return OutflowAdditionalDataHandler(lb_method.stencil, boundary_obj, target=target, field_name=field_name,
pdfs_data_type=pdfs_data_type, zeroth_timestep=zeroth_timestep)
elif isinstance(boundary_obj, NoSlipLinearBouzidi):
return NoSlipLinearBouzidiAdditionalDataHandler(lb_method.stencil, boundary_obj)
elif isinstance(boundary_obj, QuadraticBounceBack):
return QuadraticBounceBackAdditionalDataHandler(lb_method.stencil, boundary_obj)
else:
raise ValueError(f"No default AdditionalDataHandler available for boundary of type {boundary_obj.__class__}")
......@@ -107,7 +124,7 @@ class UBBAdditionalDataHandler(AdditionalDataHandler):
@property
def initialiser_list(self):
return "elementInitaliser(velocityCallback),"
return "elementInitialiser(velocityCallback),"
@property
def additional_arguments_for_fill_function(self):
......@@ -117,23 +134,117 @@ class UBBAdditionalDataHandler(AdditionalDataHandler):
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];"]
def data_initialisation(self, *_):
init_list = ["Vector3<real_t> InitialisationAdditionalData = elementInitialiser(Cell(it.x(), it.y(), it.z()), "
"blocks, *block);", "element.vel_0 = InitialisationAdditionalData[0];",
"element.vel_1 = InitialisationAdditionalData[1];"]
if self._dim == 3:
init_list.append("element.vel_2 = InitialisatonAdditionalData[2];")
init_list.append("element.vel_2 = InitialisationAdditionalData[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; "
"elementInitialiser; "
class NoSlipLinearBouzidiAdditionalDataHandler(AdditionalDataHandler):
def __init__(self, stencil, boundary_object):
assert isinstance(boundary_object, NoSlipLinearBouzidi)
self._dtype = BasicType(boundary_object.data_type).c_name
self._blocks = "const shared_ptr<StructuredBlockForest>&, IBlock&)>"
super(NoSlipLinearBouzidiAdditionalDataHandler, self).__init__(stencil=stencil)
@property
def constructor_argument_name(self):
return "wallDistanceBouzidi"
@property
def constructor_arguments(self):
return f", std::function<{self._dtype}(const Cell &, const Cell &, {self._blocks}&" \
f"{self.constructor_argument_name} "
@property
def initialiser_list(self):
return f"elementInitialiser({self.constructor_argument_name}),"
@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):
cx = self._walberla_stencil[direction][0]
cy = self._walberla_stencil[direction][1]
cz = self._walberla_stencil[direction][2]
fluid_cell = "Cell(it.x(), it.y(), it.z())"
boundary_cell = f"Cell(it.x() + {cx}, it.y() + {cy}, it.z() + {cz})"
check_str = interpolation_bc_check_template.format(cx=-cx, cy=-cy, cz=-cz)
init_element = f"elementInitialiser({fluid_cell}, {boundary_cell}, blocks, *block)"
init_list = [f"const {self._dtype} q = (({self._dtype}) {init_element});",
"element.q = q;",
check_str]
return "\n".join(init_list)
@property
def additional_member_variable(self):
return f"std::function<{self._dtype}(const Cell &, const Cell &, {self._blocks} elementInitialiser; "
class QuadraticBounceBackAdditionalDataHandler(AdditionalDataHandler):
def __init__(self, stencil, boundary_object):
assert isinstance(boundary_object, QuadraticBounceBack)
self._dtype = BasicType(boundary_object.data_type).c_name
self._blocks = "const shared_ptr<StructuredBlockForest>&, IBlock&)>"
super(QuadraticBounceBackAdditionalDataHandler, self).__init__(stencil=stencil)
@property
def constructor_argument_name(self):
return "wallDistanceQuadraticBB"
@property
def constructor_arguments(self):
return f", std::function<{self._dtype}(const Cell &, const Cell &, {self._blocks}&" \
f"{self.constructor_argument_name} "
@property
def initialiser_list(self):
return f"elementInitialiser({self.constructor_argument_name}),"
@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):
cx = self._walberla_stencil[direction][0]
cy = self._walberla_stencil[direction][1]
cz = self._walberla_stencil[direction][2]
fluid_cell = "Cell(it.x(), it.y(), it.z())"
boundary_cell = f"Cell(it.x() + {cx}, it.y() + {cy}, it.z() + {cz})"
init_element = f"elementInitialiser({fluid_cell}, {boundary_cell}, blocks, *block)"
init_list = [f"const {self._dtype} q = (({self._dtype}) {init_element});", "element.q = q;"]
return "\n".join(init_list)
@property
def additional_member_variable(self):
return f"std::function<{self._dtype}(const Cell &, const Cell &, {self._blocks} elementInitialiser; "
class OutflowAdditionalDataHandler(AdditionalDataHandler):
def __init__(self, stencil, boundary_object, target=Target.CPU, field_name='pdfs', pdfs_data_type=None, zeroth_timestep=None):
def __init__(self, stencil, boundary_object, target=Target.CPU, field_name='pdfs',
pdfs_data_type=None, zeroth_timestep=None):
assert isinstance(boundary_object, ExtrapolationOutflow)
self._stencil = boundary_object.stencil
self._lb_method = boundary_object.lb_method
......
......@@ -113,7 +113,8 @@ def __generate_alternating_lbm_boundary(generation_context,
**create_kernel_params):
if boundary_object.additional_data and additional_data_handler is None:
target = create_kernel_params.get('target', Target.CPU)
additional_data_handler = default_additional_data_handler(boundary_object, lb_method, field_name, target=target)
additional_data_handler = default_additional_data_handler(boundary_object, lb_method, field_name,
target=target, pdfs_data_type=field_data_type)
timestep_param_name = 'timestep'
timestep_param_dtype = np.uint8
......
......@@ -6,8 +6,11 @@
waLBerla_link_files_to_builddir( "*.prm" )
waLBerla_link_files_to_builddir( "*.py" )
if( WALBERLA_BUILD_WITH_CODEGEN )
waLBerla_generate_target_from_python(NAME ExampleGenerated
FILE Example.py
CODEGEN_CFG example_codegen
OUT_FILES LBMStorageSpecification.h LBMStorageSpecification.cpp
LBMSweepCollection.h LBMSweepCollection.cpp
NoSlip.h NoSlip.cpp
......@@ -16,6 +19,24 @@ waLBerla_generate_target_from_python(NAME ExampleGenerated
Example_InfoHeader.h)
waLBerla_compile_test( FILES Example.cpp DEPENDS ExampleGenerated blockforest field lbm_generated timeloop )
waLBerla_generate_target_from_python(NAME InterpolationNoSlipGenerated
FILE InterpolationNoSlip.py
CODEGEN_CFG interpolation_no_slip_codegen
OUT_FILES InterpolationNoSlipStorageSpecification.h InterpolationNoSlipStorageSpecification.cpp
InterpolationNoSlipSweepCollection.h InterpolationNoSlipSweepCollection.cpp
NoSlip.h NoSlip.cpp
NoSlipBouzidi.h NoSlipBouzidi.cpp
NoSlipQuadraticBB.h NoSlipQuadraticBB.cpp
UBB.h UBB.cpp
InterpolationNoSlipBoundaryCollection.h
InterpolationNoSlipHeader.h)
waLBerla_compile_test( FILES InterpolationNoSlip.cpp DEPENDS InterpolationNoSlipGenerated blockforest core field geometry lbm_generated timeloop )
# waLBerla_execute_test( NAME InterpolationNoSlip1 COMMAND $<TARGET_FILE:InterpolationNoSlip> ${CMAKE_CURRENT_SOURCE_DIR}/InterpolationNoSlip.prm -Parameters.distanceWall=0.1 )
# waLBerla_execute_test( NAME InterpolationNoSlip2 COMMAND $<TARGET_FILE:InterpolationNoSlip> ${CMAKE_CURRENT_SOURCE_DIR}/InterpolationNoSlip.prm -Parameters.distanceWall=0.5 )
waLBerla_execute_test( NAME InterpolationNoSlip3 COMMAND $<TARGET_FILE:InterpolationNoSlip> ${CMAKE_CURRENT_SOURCE_DIR}/InterpolationNoSlip.prm )
endif()
waLBerla_generate_target_from_python(NAME FreeSlipRefinementGenerated
FILE FreeSlipRefinement.py
CODEGEN_CFG free_slip_refinement_codegen
......
//======================================================================================================================
//
// 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 InterpolatedNoSlip.cpp
//! \author Markus Holzer <markus.holzer@fau.de>
//! \brief Couette flow driven by a UBB BC in Nord and wall boundary in the South. Remaining directions are periodic
//! If Interpolation BC are used the distance of the southern wall can be controlled. The velocity in the
//! first fluid cell is checked and compared with the velocity obtained from a default NoSlip BC.
//! Depending on the set distance for the interpolation BCs the velocity is expected to be higher or lower
//
//======================================================================================================================
#include "blockforest/Initialization.h"
#include "blockforest/communication/UniformBufferedScheme.h"
#include "core/DataTypes.h"
#include "core/Environment.h"
#include "core/debug/TestSubsystem.h"
#include "core/logging/Initialization.h"
#include "core/math/Vector3.h"
#include "core/timing/RemainingTimeLogger.h"
#include "field/AddToStorage.h"
#include "field/FlagField.h"
#include "field/GhostLayerField.h"
#include "field/vtk/VTKWriter.h"
#include "geometry/InitBoundaryHandling.h"
#include "timeloop/SweepTimeloop.h"
#include "lbm_generated/communication/UniformGeneratedPdfPackInfo.h"
#include "lbm_generated/field/AddToStorage.h"
#include "lbm_generated/field/PdfField.h"
// include the generated header file. It includes all generated classes
#include "InterpolationNoSlipHeader.h"
using namespace walberla;
using namespace std::placeholders;
using StorageSpecification_T = lbm::InterpolationNoSlipStorageSpecification;
using Stencil_T = StorageSpecification_T::Stencil;
using CommunicationStencil_T = StorageSpecification_T::CommunicationStencil;
using PdfField_T = lbm_generated::PdfField< StorageSpecification_T >;
using PackInfo_T = lbm_generated::UniformGeneratedPdfPackInfo< PdfField_T >;
using SweepCollection_T = lbm::InterpolationNoSlipSweepCollection;
using VectorField_T = GhostLayerField< real_t, StorageSpecification_T::Stencil::D >;
using ScalarField_T = GhostLayerField< real_t, 1 >;
using flag_t = walberla::uint8_t;
using FlagField_T = FlagField< flag_t >;
using BoundaryCollection_T = lbm::InterpolationNoSlipBoundaryCollection< FlagField_T >;
using blockforest::communication::UniformBufferedScheme;
class wallDistance
{
public:
wallDistance(const real_t q) : q_(q) {}
real_t operator()(const Cell& fluidCell, const Cell& boundaryCell, const shared_ptr< StructuredBlockForest >& SbF,
IBlock& block) const;
private:
const real_t q_;
}; // class wallDistance
real_t wallDistance::operator()(const Cell& /*fluidCell*/, const Cell& /*boundaryCell*/,
const shared_ptr< StructuredBlockForest >& /*SbF*/, IBlock& /*block*/) const
{
return q_;
}
int main(int argc, char** argv)
{
debug::enterTestMode();
walberla::Environment walberlaEnv(argc, argv);
logging::configureLogging(walberlaEnv.config());
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 real_t distanceWall = parameters.getParameter< real_t >("distanceWall", real_c(0.5));
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10)) + uint_c(1);
WALBERLA_LOG_DEVEL_VAR(distanceWall)
auto remainingTimeLoggerFrequency =
parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(3.0)); // in seconds
const StorageSpecification_T StorageSpec = StorageSpecification_T();
BlockDataID pdfFieldId = lbm_generated::addPdfFieldToStorage(blocks, "pdf field", StorageSpec, uint_c(1));
BlockDataID velFieldId = field::addToStorage< VectorField_T >(blocks, "Velocity", real_c(0.0), field::fzyx);
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field", uint_c(1));
SweepCollection_T sweepCollection(blocks, pdfFieldId, velFieldId, omega);
for (auto& block : *blocks)
{
sweepCollection.initialise(&block);
}
const FlagUID fluidFlagUID("Fluid");
auto boundariesConfig = walberlaEnv.config()->getBlock("Boundaries");
geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldId, boundariesConfig);
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID);
const wallDistance wallDistanceCallback{ distanceWall };
std::function< real_t(const Cell&, const Cell&, const shared_ptr< StructuredBlockForest >&, IBlock&) >
wallDistanceFunctor = wallDistanceCallback;
// For the BoundaryCollection a funcotr to calculate the wall distance for the Bouzidi NoSlip and for the QuadraticBB
// have to be provided. In this test case we use the same function to calculate the wall distance
BoundaryCollection_T boundaryCollection(blocks, flagFieldId, pdfFieldId, fluidFlagUID, omega, wallDistanceFunctor,
wallDistanceFunctor);
auto packInfo = std::make_shared< lbm_generated::UniformGeneratedPdfPackInfo< PdfField_T > >(pdfFieldId);
UniformBufferedScheme< Stencil_T > communication(blocks);
communication.addPackInfo(packInfo);
SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
timeloop.add() << BeforeFunction(communication, "communication")
<< Sweep(boundaryCollection.getSweep(BoundaryCollection_T::ALL), "Boundary Conditions");
timeloop.add() << Sweep(sweepCollection.streamCollide(SweepCollection_T::ALL), "LBM StreamCollide");
uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0);
if (vtkWriteFrequency > 0)
{
auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "InterpolationNoSlipVTK", vtkWriteFrequency, 0, false,
"vtk_out", "simulation_step", false, true, true, false, 0);
auto velWriter = make_shared< field::VTKWriter< VectorField_T > >(velFieldId, "velocity");
vtkOutput->addBeforeFunction([&]() {
for (auto& block : *blocks)
{
sweepCollection.calculateMacroscopicParameters(&block);
}
});
vtkOutput->addCellDataWriter(velWriter);
timeloop.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
}
if (remainingTimeLoggerFrequency > 0)
{
// log remaining time
timeloop.addFuncAfterTimeStep(
timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
"remaining time logger");
}
timeloop.run();
// This is the velocity at the wall, when a NoSlip BC is used. This is similar to using an interpolation BC with a
// wall distance of 0.5. This value can be obtained by either setting distanceWall to 0.5 in the Parameter file or
// specifying the NoSlip BC at the southern boundary
const real_t defaultNoSlipVelocity = real_c(0.0002);
for (auto& block : *blocks)
{
sweepCollection.calculateMacroscopicParameters(&block);
auto velField = block.getData< VectorField_T >(velFieldId);
auto velAtWall = velField->get(cell_idx_c(cellsPerBlock[0] / 2), 0, cell_idx_c(cellsPerBlock[2] / 2), 0);
// WALBERLA_LOG_DEVEL_VAR(velAtWall)
if (distanceWall > 0.49 && distanceWall < 0.51) { WALBERLA_CHECK_FLOAT_EQUAL(velAtWall, defaultNoSlipVelocity) }
else if (distanceWall < 0.49) { WALBERLA_CHECK_GREATER(defaultNoSlipVelocity, velAtWall) }
else if (distanceWall > 0.51) { WALBERLA_CHECK_LESS(defaultNoSlipVelocity, velAtWall) }
}
return EXIT_SUCCESS;
}
Parameters
{
omega 1.1;
timesteps 5000;
distanceWall 0.9;
remainingTimeLoggerFrequency 0; // in seconds
vtkWriteFrequency 0;
}
DomainSetup
{
blocks < 1, 1, 1 >;
cellsPerBlock < 50, 25, 25 >;
periodic < 1, 0, 1 >;
}
Boundaries
{
// Border { direction S; walldistance -1; flag NoSlip; }
// Border { direction S; walldistance -1; flag NoSlipBouzidi; }
Border { direction S; walldistance -1; flag NoSlipQuadraticBB; }
Border { direction N; walldistance -1; flag UBB; }
}
Logging
{
logLevel info; // info progress detail tracing
}
import sympy as sp
from pystencils import Target
from pystencils import fields
from lbmpy.advanced_streaming.utility import get_timesteps
from lbmpy.boundaries import NoSlip, NoSlipLinearBouzidi, QuadraticBounceBack, UBB
from lbmpy.creationfunctions import create_lb_method, create_lb_collision_rule
from lbmpy import LBMConfig, LBMOptimisation, Stencil, Method, LBStencil
from pystencils_walberla import CodeGeneration, generate_info_header
from lbmpy_walberla import generate_lbm_package, lbm_boundary_generator
import warnings
warnings.filterwarnings("ignore")
with CodeGeneration() as ctx:
target = Target.CPU # Target.GPU if ctx.cuda else Target.CPU
data_type = "float64" if ctx.double_accuracy else "float32"
streaming_pattern = 'pull'
timesteps = get_timesteps(streaming_pattern)
omega = sp.symbols("omega")
stencil = LBStencil(Stencil.D3Q27)
pdfs, vel_field = fields(f"pdfs({stencil.Q}), velocity({stencil.D}): {data_type}[{stencil.D}D]",
layout='fzyx')
macroscopic_fields = {'velocity': vel_field}
lbm_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=omega,
streaming_pattern=streaming_pattern)
lbm_opt = LBMOptimisation(cse_global=False, field_layout='fzyx')
method = create_lb_method(lbm_config=lbm_config)
collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
no_slip = lbm_boundary_generator(class_name='NoSlip', flag_uid='NoSlip',
boundary_object=NoSlip())
no_slip_bouzidi = lbm_boundary_generator(class_name='NoSlipBouzidi', flag_uid='NoSlipBouzidi',
boundary_object=NoSlipLinearBouzidi(data_type=data_type))
no_slip_quadraticbb = lbm_boundary_generator(class_name='NoSlipQuadraticBB', flag_uid='NoSlipQuadraticBB',
boundary_object=QuadraticBounceBack(omega, data_type=data_type))
ubb = lbm_boundary_generator(class_name='UBB', flag_uid='UBB',
boundary_object=UBB([0.01, 0, 0], data_type=data_type))
generate_lbm_package(ctx, name="InterpolationNoSlip",
collision_rule=collision_rule,
lbm_config=lbm_config, lbm_optimisation=lbm_opt,
nonuniform=True, boundaries=[no_slip, no_slip_bouzidi, no_slip_quadraticbb, ubb],
macroscopic_fields=macroscopic_fields, data_type=data_type)
generate_info_header(ctx, 'InterpolationNoSlipHeader')
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