Skip to content
Snippets Groups Projects
Commit 2b47f2a7 authored by Philipp Suffa's avatar Philipp Suffa
Browse files

Merge branch 'InterpolcationBCIntegration' into 'master'

Integration for Interpolation BCs

See merge request walberla/walberla!652
parents b346f405 d514fba4
Branches
No related tags found
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% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment