Commit 0e4b7712 authored by Frederik Hennig's avatar Frederik Hennig Committed by Markus Holzer
Browse files

Code Generation for Lattice Boltzmann In-Place Streaming Patterns

parent d884044c
......@@ -21,6 +21,7 @@ if ( WALBERLA_BUILD_WITH_PYTHON )
add_subdirectory( FieldCommunication )
if ( WALBERLA_BUILD_WITH_CODEGEN )
add_subdirectory( FlowAroundSphereCodeGen )
add_subdirectory( UniformGridGenerated )
add_subdirectory( PhaseFieldAllenCahn )
endif()
......
waLBerla_link_files_to_builddir( "*.py" )
if (WALBERLA_BUILD_WITH_CUDA)
waLBerla_generate_target_from_python(NAME FlowAroundSphereGenerated
FILE FlowAroundSphereCodeGen.py
OUT_FILES FlowAroundSphereCodeGen_LbSweep.cu FlowAroundSphereCodeGen_LbSweep.h
FlowAroundSphereCodeGen_MacroSetter.cu FlowAroundSphereCodeGen_MacroSetter.h
FlowAroundSphereCodeGen_UBB.cu FlowAroundSphereCodeGen_UBB.h
FlowAroundSphereCodeGen_NoSlip.cu FlowAroundSphereCodeGen_NoSlip.h
FlowAroundSphereCodeGen_Outflow.cu FlowAroundSphereCodeGen_Outflow.h
FlowAroundSphereCodeGen_PackInfoEven.cu FlowAroundSphereCodeGen_PackInfoEven.h
FlowAroundSphereCodeGen_PackInfoOdd.cu FlowAroundSphereCodeGen_PackInfoOdd.h
FlowAroundSphereCodeGen_InfoHeader.h)
waLBerla_add_executable( NAME FlowAroundSphereCodeGen FILE FlowAroundSphereCodeGen.cpp
DEPENDS blockforest boundary core domain_decomposition field geometry timeloop vtk FlowAroundSphereGenerated)
set_target_properties( FlowAroundSphereCodeGen PROPERTIES CXX_VISIBILITY_PRESET hidden)
else ()
waLBerla_generate_target_from_python(NAME FlowAroundSphereGenerated
FILE FlowAroundSphereCodeGen.py
OUT_FILES FlowAroundSphereCodeGen_LbSweep.cpp FlowAroundSphereCodeGen_LbSweep.h
FlowAroundSphereCodeGen_MacroSetter.cpp FlowAroundSphereCodeGen_MacroSetter.h
FlowAroundSphereCodeGen_UBB.cpp FlowAroundSphereCodeGen_UBB.h
FlowAroundSphereCodeGen_NoSlip.cpp FlowAroundSphereCodeGen_NoSlip.h
FlowAroundSphereCodeGen_Outflow.cpp FlowAroundSphereCodeGen_Outflow.h
FlowAroundSphereCodeGen_PackInfoEven.cpp FlowAroundSphereCodeGen_PackInfoEven.h
FlowAroundSphereCodeGen_PackInfoOdd.cpp FlowAroundSphereCodeGen_PackInfoOdd.h
FlowAroundSphereCodeGen_InfoHeader.h)
waLBerla_add_executable( NAME FlowAroundSphereCodeGen FILE FlowAroundSphereCodeGen.cpp
DEPENDS blockforest boundary core domain_decomposition field geometry timeloop vtk FlowAroundSphereGenerated)
set_target_properties( FlowAroundSphereCodeGen PROPERTIES CXX_VISIBILITY_PRESET hidden)
endif()
\ No newline at end of file
//======================================================================================================================
//
// 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 FlowAroundSphereCodeGen.cpp
//! \author Frederik Hennig <frederik.hennig@fau.de>
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#include "blockforest/all.h"
#include "core/all.h"
#include "domain_decomposition/all.h"
#include "field/all.h"
#include "geometry/all.h"
#include "lbm/inplace_streaming/TimestepTracker.h"
#include "lbm/vtk/QCriterion.h"
#include "python_coupling/CreateConfig.h"
#include "python_coupling/PythonCallback.h"
#include "timeloop/all.h"
#if defined(WALBERLA_BUILD_WITH_CUDA)
# include "cuda/AddGPUFieldToStorage.h"
# include "cuda/DeviceSelectMPI.h"
# include "cuda/HostFieldAllocator.h"
# include "cuda/NVTX.h"
# include "cuda/ParallelStreams.h"
# include "cuda/communication/GPUPackInfo.h"
# include "cuda/communication/UniformGPUScheme.h"
#endif
// CodeGen includes
#include "FlowAroundSphereCodeGen_InfoHeader.h"
namespace walberla
{
typedef lbm::FlowAroundSphereCodeGen_PackInfoEven PackInfoEven_T;
typedef lbm::FlowAroundSphereCodeGen_PackInfoOdd PackInfoOdd_T;
typedef walberla::uint8_t flag_t;
typedef FlagField< flag_t > FlagField_T;
#if defined(WALBERLA_BUILD_WITH_CUDA)
typedef cuda::GPUField< real_t > GPUField;
#endif
using namespace std::placeholders;
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 > >());
};
auto VelocityCallback = [](const Cell& pos, const shared_ptr< StructuredBlockForest >& SbF, IBlock& block,
real_t inflow_velocity) {
Cell globalCell;
CellInterval domain = SbF->getDomainCellBB();
real_t h_y = domain.yMax() - domain.yMin();
real_t h_z = domain.zMax() - domain.zMin();
SbF->transformBlockLocalToGlobalCell(globalCell, block, pos);
real_t y1 = globalCell[1] - (h_y / 2.0 + 0.5);
real_t z1 = globalCell[2] - (h_z / 2.0 + 0.5);
real_t u = (inflow_velocity * 16) / (h_y * h_y * h_z * h_z) * (h_y / 2.0 - y1) * (h_y / 2 + y1) * (h_z / 2 - z1) *
(h_z / 2 + z1);
Vector3< real_t > result(u, 0.0, 0.0);
return result;
};
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_;
};
class Filter
{
public:
explicit Filter(Vector3< uint_t > numberOfCells) : numberOfCells_(numberOfCells) {}
void operator()(const IBlock& /*block*/) {}
bool operator()(const cell_idx_t x, const cell_idx_t y, const cell_idx_t z) const
{
return x >= -1 && x <= cell_idx_t(numberOfCells_[0]) && y >= -1 && y <= cell_idx_t(numberOfCells_[1]) &&
z >= -1 && z <= cell_idx_t(numberOfCells_[2]);
}
private:
Vector3< uint_t > numberOfCells_;
};
using FluidFilter_T = Filter;
int main(int argc, char** argv)
{
walberla::Environment walberlaEnv(argc, argv);
#if defined(WALBERLA_BUILD_WITH_CUDA)
cuda::selectDeviceBasedOnMpiRank();
#endif
for (auto cfg = python_coupling::configBegin(argc, argv); cfg != python_coupling::configEnd(); ++cfg)
{
WALBERLA_MPI_WORLD_BARRIER()
auto config = *cfg;
logging::configureLogging(config);
auto blocks = blockforest::createUniformBlockGridFromConfig(config);
// read parameters
Vector3< uint_t > cellsPerBlock =
config->getBlock("DomainSetup").getParameter< Vector3< uint_t > >("cellsPerBlock");
auto parameters = config->getOneBlock("Parameters");
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
const real_t omega = parameters.getParameter< real_t >("omega", real_t(1.9));
const real_t u_max = parameters.getParameter< real_t >("u_max", real_t(0.05));
const real_t reynolds_number = parameters.getParameter< real_t >("reynolds_number", real_t(1000));
const uint_t diameter_sphere = parameters.getParameter< uint_t >("diameter_sphere", uint_t(5));
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(0), field::fzyx);
#if defined(WALBERLA_BUILD_WITH_CUDA)
BlockDataID pdfFieldIDGPU = cuda::addGPUFieldToStorage< PdfField_T >(blocks, pdfFieldID, "PDFs on GPU", true);
BlockDataID velFieldIDGPU =
cuda::addGPUFieldToStorage< VelocityField_T >(blocks, velFieldID, "velocity on GPU", true);
BlockDataID densityFieldIDGPU =
cuda::addGPUFieldToStorage< ScalarField_T >(blocks, densityFieldID, "density on GPU", true);
#endif
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
// initialise all PDFs
#if defined(WALBERLA_BUILD_WITH_CUDA)
pystencils::FlowAroundSphereCodeGen_MacroSetter setterSweep(pdfFieldIDGPU, velFieldIDGPU);
for (auto& block : *blocks)
setterSweep(&block);
cuda::fieldCpy< PdfField_T, GPUField >(blocks, pdfFieldID, pdfFieldIDGPU);
#else
pystencils::FlowAroundSphereCodeGen_MacroSetter setterSweep(pdfFieldID, velFieldID);
for (auto& block : *blocks)
setterSweep(&block);
#endif
// Create communication
#if defined(WALBERLA_BUILD_WITH_CUDA)
// This way of using alternating pack infos is temporary and will soon be replaced
// by something more straight-forward
cuda::communication::UniformGPUScheme< Stencil_T > comEven(blocks, false);
comEven.addPackInfo(make_shared< PackInfoEven_T >(pdfFieldIDGPU));
auto evenComm = std::function< void() >([&]() { comEven.communicate(nullptr); });
cuda::communication::UniformGPUScheme< Stencil_T > comODD(blocks, false);
comODD.addPackInfo(make_shared< PackInfoOdd_T >(pdfFieldIDGPU));
auto oddComm = std::function< void() >([&]() { comODD.communicate(nullptr); });
#else
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));
#endif
// create and initialize boundary handling
const FlagUID fluidFlagUID("Fluid");
auto boundariesConfig = config->getOneBlock("Boundaries");
std::function< Vector3< real_t >(const Cell&, const shared_ptr< StructuredBlockForest >&, IBlock&) >
velocity_initialisation = std::bind(VelocityCallback, _1, _2, _3, u_max);
#if defined(WALBERLA_BUILD_WITH_CUDA)
lbm::FlowAroundSphereCodeGen_UBB ubb(blocks, pdfFieldIDGPU, velocity_initialisation);
lbm::FlowAroundSphereCodeGen_NoSlip noSlip(blocks, pdfFieldIDGPU);
lbm::FlowAroundSphereCodeGen_Outflow outflow(blocks, pdfFieldIDGPU, pdfFieldID);
lbm::FlowAroundSphereCodeGen_LbSweep lbSweep(densityFieldIDGPU, pdfFieldIDGPU, velFieldIDGPU, omega);
#else
lbm::FlowAroundSphereCodeGen_UBB ubb(blocks, pdfFieldID, velocity_initialisation);
lbm::FlowAroundSphereCodeGen_NoSlip noSlip(blocks, pdfFieldID);
lbm::FlowAroundSphereCodeGen_Outflow outflow(blocks, pdfFieldID);
lbm::FlowAroundSphereCodeGen_LbSweep lbSweep(densityFieldID, pdfFieldID, velFieldID, omega);
#endif
geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldId, boundariesConfig);
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID);
ubb.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("UBB"), fluidFlagUID);
noSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("NoSlip"), fluidFlagUID);
outflow.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("Outflow"), fluidFlagUID);
// create time loop
SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
// Timestep Tracking and Sweeps
auto tracker = make_shared< lbm::TimestepTracker >(0);
AlternatingBeforeFunction communication(evenComm, oddComm, tracker);
// add LBM sweep and communication to time loop
timeloop.add() << BeforeFunction(communication, "communication")
<< Sweep(noSlip.getSweep(tracker), "noSlip boundary");
timeloop.add() << Sweep(outflow.getSweep(tracker), "outflow boundary");
timeloop.add() << Sweep(ubb.getSweep(tracker), "ubb boundary");
timeloop.add() << BeforeFunction(tracker->getAdvancementFunction(), "Timestep Advancement")
<< Sweep(lbSweep.getSweep(tracker), "LB update rule");
// LBM stability check
timeloop.addFuncAfterTimeStep(makeSharedFunctor(field::makeStabilityChecker< PdfField_T, FlagField_T >(
config, blocks, pdfFieldID, flagFieldId, fluidFlagUID)),
"LBM stability check");
// log remaining time
timeloop.addFuncAfterTimeStep(
timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
"remaining time logger");
// add VTK output to time loop
uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0);
if (vtkWriteFrequency > 0)
{
auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false, "vtk_out",
"simulation_step", false, true, true, false, 0);
#if defined(WALBERLA_BUILD_WITH_CUDA)
vtkOutput->addBeforeFunction([&]() {
cuda::fieldCpy< VelocityField_T, GPUField >(blocks, velFieldID, velFieldIDGPU);
cuda::fieldCpy< ScalarField_T, GPUField >(blocks, densityFieldID, densityFieldIDGPU);
});
#endif
auto velWriter = make_shared< field::VTKWriter< VelocityField_T > >(velFieldID, "velocity");
auto densityWriter = make_shared< field::VTKWriter< ScalarField_T > >(densityFieldID, "density");
FluidFilter_T filter(cellsPerBlock);
auto QCriterionWriter = make_shared< lbm::QCriterionVTKWriter< VelocityField_T, FluidFilter_T > >(
blocks, filter, velFieldID, "Q-Criterion");
vtkOutput->addCellDataWriter(velWriter);
vtkOutput->addCellDataWriter(densityWriter);
vtkOutput->addCellDataWriter(QCriterionWriter);
timeloop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
}
WcTimer simTimer;
WALBERLA_LOG_INFO_ON_ROOT("Simulating flow around sphere:"
"\n timesteps: "
<< timesteps << "\n reynolds number: " << reynolds_number
<< "\n relaxation rate: " << omega << "\n maximum inflow velocity: " << u_max
<< "\n diameter_sphere: " << diameter_sphere)
simTimer.start();
timeloop.run();
simTimer.end();
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
auto time = simTimer.last();
auto nrOfCells = real_c(cellsPerBlock[0] * cellsPerBlock[1] * cellsPerBlock[2]);
auto mlupsPerProcess = nrOfCells * real_c(timesteps) / time * 1e-6;
WALBERLA_LOG_RESULT_ON_ROOT("MLUPS per process " << mlupsPerProcess)
WALBERLA_LOG_RESULT_ON_ROOT("Time per time step " << time / real_c(timesteps))
}
return EXIT_SUCCESS;
}
} // namespace walberla
int main(int argc, char** argv) { walberla::main(argc, argv); }
from pystencils.field import fields
from lbmpy.advanced_streaming.utility import get_timesteps, Timestep
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
from lbmpy.stencils import get_stencil
from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_method, create_lb_update_rule
from lbmpy.boundaries import NoSlip, UBB, ExtrapolationOutflow
from pystencils_walberla import CodeGeneration, generate_sweep, generate_info_header
from lbmpy_walberla.additional_data_handler import UBBAdditionalDataHandler, OutflowAdditionalDataHandler
from lbmpy_walberla import generate_boundary, generate_lb_pack_info
from lbmpy_walberla import generate_alternating_lbm_sweep, generate_alternating_lbm_boundary
import sympy as sp
stencil = get_stencil("D3Q27")
q = len(stencil)
dim = len(stencil[0])
streaming_pattern = 'esotwist'
timesteps = get_timesteps(streaming_pattern)
pdfs, velocity_field, density_field = fields(f"pdfs({q}), velocity({dim}), density(1) : double[{dim}D]", layout='fzyx')
omega = sp.Symbol("omega")
u_max = sp.Symbol("u_max")
output = {
'density': density_field,
'velocity': velocity_field
}
opt = {'symbolic_field': pdfs,
'cse_global': False,
'cse_pdfs': False}
method_params = {'method': 'cumulant',
'stencil': stencil,
'relaxation_rate': omega,
'galilean_correction': True,
'field_name': 'pdfs',
'streaming_pattern': streaming_pattern,
'output': output,
'optimization': opt}
collision_rule = create_lb_collision_rule(**method_params)
lb_method = collision_rule.method
# getter & setter
setter_assignments = macroscopic_values_setter(lb_method, velocity=velocity_field.center_vector,
pdfs=pdfs, density=1,
streaming_pattern=streaming_pattern,
previous_timestep=timesteps[0])
# opt = {'instruction_set': 'sse', 'assume_aligned': True, 'nontemporal': False, 'assume_inner_stride_one': True}
stencil_typedefs = {'Stencil_T': stencil}
field_typedefs = {'PdfField_T': pdfs,
'VelocityField_T': velocity_field,
'ScalarField_T': density_field}
with CodeGeneration() as ctx:
if ctx.cuda:
target = 'gpu'
else:
target = 'cpu'
opt['target'] = target
# sweeps
generate_alternating_lbm_sweep(ctx, 'FlowAroundSphereCodeGen_LbSweep',
collision_rule, streaming_pattern, optimization=opt)
generate_sweep(ctx, 'FlowAroundSphereCodeGen_MacroSetter', setter_assignments, target=target)
# boundaries
ubb = UBB(lambda *args: None, dim=dim)
ubb_data_handler = UBBAdditionalDataHandler(stencil, ubb)
outflow = ExtrapolationOutflow(stencil[4], lb_method, streaming_pattern=streaming_pattern)
outflow_data_handler = OutflowAdditionalDataHandler(stencil, outflow, target=target)
generate_alternating_lbm_boundary(ctx, 'FlowAroundSphereCodeGen_UBB', ubb, lb_method,
target=target, streaming_pattern=streaming_pattern,
additional_data_handler=ubb_data_handler)
generate_alternating_lbm_boundary(ctx, 'FlowAroundSphereCodeGen_NoSlip', NoSlip(), lb_method,
target=target, streaming_pattern=streaming_pattern)
generate_alternating_lbm_boundary(ctx, 'FlowAroundSphereCodeGen_Outflow', outflow, lb_method,
target=target, streaming_pattern=streaming_pattern,
additional_data_handler=outflow_data_handler)
# communication
generate_lb_pack_info(ctx, 'FlowAroundSphereCodeGen_PackInfo', stencil, pdfs,
streaming_pattern=streaming_pattern, always_generate_separate_classes=True, target=target)
# Info header containing correct template definitions for stencil and field
generate_info_header(ctx, 'FlowAroundSphereCodeGen_InfoHeader',
stencil_typedefs=stencil_typedefs, field_typedefs=field_typedefs)
import waLBerla as wlb
from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
class Scenario:
def __init__(self):
self.timesteps = 101
self.vtkWriteFrequency = 2000
self.cells = (384, 128, 128)
self.blocks = (1, 1, 1)
self.periodic = (0, 0, 0)
self.diameter_sphere = min(self.cells) // 2
self.u_max = 0.1
self.reynolds_number = 1000000
kinematic_viscosity = (self.diameter_sphere * self.u_max) / self.reynolds_number
self.omega = relaxation_rate_from_lattice_viscosity(kinematic_viscosity)
self.total_cells = (self.cells[0] * self.blocks[0],
self.cells[1] * self.blocks[1],
self.cells[2] * self.blocks[2])
@wlb.member_callback
def config(self):
return {
'DomainSetup': {
'blocks': self.blocks,
'cellsPerBlock': self.cells,
'periodic': self.periodic,
'oneBlockPerProcess': True
},
'Parameters': {
'timesteps': self.timesteps,
'vtkWriteFrequency': self.vtkWriteFrequency,
'omega': self.omega,
'u_max': self.u_max,
'reynolds_number': self.reynolds_number,
'diameter_sphere': self.diameter_sphere
},
'Boundaries': {
'Border': [
{'direction': 'N', 'walldistance': -1, 'flag': 'NoSlip'},
{'direction': 'S', 'walldistance': -1, 'flag': 'NoSlip'},
{'direction': 'W', 'walldistance': -1, 'flag': 'UBB'},
{'direction': 'E', 'walldistance': -1, 'flag': 'Outflow'},
{'direction': 'T', 'walldistance': -1, 'flag': 'NoSlip'},
{'direction': 'B', 'walldistance': -1, 'flag': 'NoSlip'},
],
'Body': [
{'shape': "sphere",
'midpoint': (int(0.40 * self.total_cells[0]), self.total_cells[1] // 2, self.total_cells[2] // 2),
'radius': self.diameter_sphere // 2,
'flag': 'NoSlip'}
]
},
}
scenarios = wlb.ScenarioManager()
scenarios.add(Scenario())
from .boundary import generate_boundary
from .boundary import generate_boundary, generate_alternating_lbm_boundary
from .walberla_lbm_generation import RefinementScaling, generate_lattice_model
from .packinfo import generate_lb_pack_info
from .alternating_sweeps import generate_alternating_lbm_sweep
__all__ = ['generate_lattice_model', 'RefinementScaling', 'generate_boundary', 'generate_lb_pack_info']
__all__ = ['generate_lattice_model', 'generate_alternating_lbm_sweep',
'RefinementScaling', 'generate_boundary', 'generate_alternating_lbm_boundary',
'generate_lb_pack_info']
from pystencils.stencil import inverse_direction
from lbmpy.advanced_streaming import AccessPdfValues, numeric_offsets, numeric_index
from lbmpy.boundaries.boundaryconditions import LbBoundary
from lbmpy.boundaries import ExtrapolationOutflow, UBB
from pystencils_walberla.additional_data_handler import AdditionalDataHandler
def default_additional_data_handler(boundary_obj: LbBoundary, lb_method, field_name, target='cpu'):
if not boundary_obj.additional_data:
return None
if 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)
else:
raise ValueError(f"No default AdditionalDataHandler available for boundary of type {boundary_obj.__class__}")
class UBBAdditionalDataHandler(AdditionalDataHandler):
def __init__(self, stencil, boundary_object):
assert isinstance(boundary_object, UBB)
......@@ -114,7 +127,7 @@ class OutflowAdditionalDataHandler(AdditionalDataHandler):
if self._dim == 2:
pos.append("0")
pos.append(str(numeric_index(pdf_accessor.accs[inv_dir])[0]))
result[f'pdf'] = ', '.join(pos)
result['pdf'] = ', '.join(pos)
pos = []
for p, o, t in zip(position, offsets, tangential_offset):
......@@ -122,6 +135,6 @@ class OutflowAdditionalDataHandler(AdditionalDataHandler):
if self._dim == 2:
pos.append("0")
pos.append(str(numeric_index(pdf_accessor.accs[inv_dir])[0]))
result[f'pdf_nd'] = ', '.join(pos)
result['pdf_nd'] = ', '.join(pos)
return result
import numpy as np
from pystencils_walberla.codegen import generate_selective_sweep, get_vectorize_instruction_set
from pystencils_walberla.kernel_selection import (
AbstractInterfaceArgumentMapping, AbstractConditionNode, KernelCallNode)
from pystencils import TypedSymbol