From a609c01191fb806376633ab63e96273a02354e5c Mon Sep 17 00:00:00 2001 From: Markus Holzer Date: Wed, 31 Jan 2024 14:48:51 +0100 Subject: [PATCH] FreeSlip BC with refinement --- .../templates/Boundary.tmpl.h | 1 + .../refinement/BasicRecursiveTimeStep.impl.h | 3 +- tests/lbm_generated/CMakeLists.txt | 12 + tests/lbm_generated/FreeSlipRefinement.cpp | 280 ++++++++++++++++++ tests/lbm_generated/FreeSlipRefinement.prm | 26 ++ tests/lbm_generated/FreeSlipRefinement.py | 49 +++ 6 files changed, 369 insertions(+), 2 deletions(-) create mode 100644 tests/lbm_generated/FreeSlipRefinement.cpp create mode 100644 tests/lbm_generated/FreeSlipRefinement.prm create mode 100644 tests/lbm_generated/FreeSlipRefinement.py diff --git a/python/pystencils_walberla/templates/Boundary.tmpl.h b/python/pystencils_walberla/templates/Boundary.tmpl.h index 75e3cd13..704a7227 100644 --- a/python/pystencils_walberla/templates/Boundary.tmpl.h +++ b/python/pystencils_walberla/templates/Boundary.tmpl.h @@ -19,6 +19,7 @@ #pragma once #include "core/DataTypes.h" +#include "core/logging/Logging.h" {% if target is equalto 'cpu' -%} #include "field/GhostLayerField.h" diff --git a/src/lbm_generated/refinement/BasicRecursiveTimeStep.impl.h b/src/lbm_generated/refinement/BasicRecursiveTimeStep.impl.h index 215f5a5c..3abb8a91 100644 --- a/src/lbm_generated/refinement/BasicRecursiveTimeStep.impl.h +++ b/src/lbm_generated/refinement/BasicRecursiveTimeStep.impl.h @@ -178,8 +178,7 @@ std::function BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, Bo template< typename PdfField_T, typename SweepCollection_T, typename BoundaryCollection_T > -void BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T >::ghostLayerPropagation( - Block * block) +void BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T >::ghostLayerPropagation(Block * block) { auto pdfField = block->getData(pdfFieldId_); diff --git a/tests/lbm_generated/CMakeLists.txt b/tests/lbm_generated/CMakeLists.txt index d7a7ef76..c221e14c 100644 --- a/tests/lbm_generated/CMakeLists.txt +++ b/tests/lbm_generated/CMakeLists.txt @@ -16,6 +16,18 @@ 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 FreeSlipRefinementGenerated + FILE FreeSlipRefinement.py + CODEGEN_CFG free_slip_refinement_codegen + OUT_FILES FreeSlipRefinementStorageSpecification.h FreeSlipRefinementStorageSpecification.cpp + FreeSlipRefinementSweepCollection.h FreeSlipRefinementSweepCollection.cpp + FreeSlip.h FreeSlip.cpp + UBB.h UBB.cpp + Outflow.h Outflow.cpp + FreeSlipRefinementBoundaryCollection.h + FreeSlipRefinementInfoHeader.h) +waLBerla_compile_test( FILES FreeSlipRefinement.cpp DEPENDS FreeSlipRefinementGenerated blockforest field lbm_generated timeloop ) + if( WALBERLA_DOUBLE_ACCURACY ) waLBerla_compile_test( FILES LDC.cpp DEPENDS blockforest field lbm_generated timeloop ) endif() diff --git a/tests/lbm_generated/FreeSlipRefinement.cpp b/tests/lbm_generated/FreeSlipRefinement.cpp new file mode 100644 index 00000000..42d91a6e --- /dev/null +++ b/tests/lbm_generated/FreeSlipRefinement.cpp @@ -0,0 +1,280 @@ +//====================================================================================================================== +// +// 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 . +// +//! \file FreeSlipRefinement.cpp +//! \author Markus Holzer +//! \brief Channel flow with inlet and outlet on West and East. The rest of the Boundaries consist of FreeSlip. The +//! Channel flow should reach the inlet velocity in the whole domain because the FreeSlip BC will not provide a +//! restriction on it. With the D3Q27 stencil this only works if the BC is also set on the first fluid node +// +//====================================================================================================================== + +#include "blockforest/Initialization.h" +#include "blockforest/communication/NonUniformBufferedScheme.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/NonuniformGeneratedPdfPackInfo.h" +#include "lbm_generated/field/AddToStorage.h" +#include "lbm_generated/field/PdfField.h" +#include "lbm_generated/refinement/BasicRecursiveTimeStep.h" + +// include the generated header file. It includes all generated classes +#include "FreeSlipRefinementInfoHeader.h" + +using namespace walberla; +using namespace std::placeholders; + +using StorageSpecification_T = lbm::FreeSlipRefinementStorageSpecification; +using Stencil_T = StorageSpecification_T::Stencil; +using CommunicationStencil_T = StorageSpecification_T::CommunicationStencil; +using PdfField_T = lbm_generated::PdfField< StorageSpecification_T >; + +using SweepCollection_T = lbm::FreeSlipRefinementSweepCollection; + +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::FreeSlipRefinementBoundaryCollection< FlagField_T >; + +using RefinementSelectionFunctor = SetupBlockForest::RefinementSelectionFunction; + +class ChannelRefinement +{ + public: + ChannelRefinement(const uint_t depth) : refinementDepth_(depth){}; + + void operator()(SetupBlockForest& forest) + { + std::vector< SetupBlock* > blocks; + forest.getBlocks(blocks); + AABB refinementAABB = AABB(forest.getDomain().xSize() / 2 - 1, forest.getDomain().yMin(), forest.getDomain().zSize() / 2 - 1, + forest.getDomain().xSize() / 2 + 1, forest.getDomain().yMin() + 1, forest.getDomain().zSize() / 2 + 1 ); + + for (auto b : blocks) + { + if (refinementAABB.intersects(b->getAABB()) && b->getLevel() < refinementDepth_) + { + b->setMarker(true); + } + } + } + + private: + const uint_t refinementDepth_; +}; + +class Channel +{ + public: + Channel(const uint_t depth) : refinementDepth_(depth), freeSlipFlagUID_("FreeSlip"), ubbFlagUID_("UBB"), outflowFlagUID_("Outflow"){}; + + RefinementSelectionFunctor refinementSelector() { return ChannelRefinement(refinementDepth_); } + void setupBoundaryFlagField(StructuredBlockForest& sbfs, const BlockDataID flagFieldID) + { + for (auto bIt = sbfs.begin(); bIt != sbfs.end(); ++bIt) + { + Block& b = dynamic_cast< Block& >(*bIt); + uint_t level = b.getLevel(); + auto flagField = b.getData< FlagField_T >(flagFieldID); + uint8_t freeSlipFlag = flagField->registerFlag(freeSlipFlagUID_); + uint8_t ubbFlag = flagField->registerFlag(ubbFlagUID_); + uint8_t outflowFlag = flagField->registerFlag(outflowFlagUID_); + for (auto cIt = flagField->beginWithGhostLayerXYZ(2); cIt != flagField->end(); ++cIt) + { + Cell localCell = cIt.cell(); + Cell globalCell(localCell); + sbfs.transformBlockLocalToGlobalCell(globalCell, b); + if (globalCell.x() < 0) { flagField->addFlag(localCell, ubbFlag); } + else if (globalCell.x() >= cell_idx_c(sbfs.getNumberOfXCells(level))) { flagField->addFlag(localCell, outflowFlag); } + else if (globalCell.y() >= cell_idx_c(sbfs.getNumberOfYCells(level))) { flagField->addFlag(localCell, freeSlipFlag); } + else if (globalCell.y() < cell_idx_c(1 << level)) {flagField->addFlag(localCell, freeSlipFlag);} + } + } + } + + private: + const std::string refinementProfile_; + const uint_t refinementDepth_; + + const FlagUID freeSlipFlagUID_; + const FlagUID ubbFlagUID_; + const FlagUID outflowFlagUID_; +}; + +static void createSetupBlockForest(SetupBlockForest& setupBfs, const Config::BlockHandle& domainSetup, Channel& setup) +{ + Vector3< real_t > domainSize = domainSetup.getParameter< Vector3< real_t > >("domainSize"); + Vector3< uint_t > rootBlocks = domainSetup.getParameter< Vector3< uint_t > >("rootBlocks"); + Vector3< bool > periodic = domainSetup.getParameter< Vector3< bool > >("periodic"); + + auto refSelection = setup.refinementSelector(); + setupBfs.addRefinementSelectionFunction(std::function< void(SetupBlockForest&) >(refSelection)); + AABB domain(real_t(0.0), real_t(0.0), real_t(0.0), domainSize[0], domainSize[1], domainSize[2]); + setupBfs.init(domain, rootBlocks[0], rootBlocks[1], rootBlocks[2], periodic[0], periodic[1], periodic[2]); + setupBfs.balanceLoad(blockforest::StaticLevelwiseCurveBalance(true), uint_c(MPIManager::instance()->numProcesses())); +} + +int main(int argc, char** argv) +{ + walberla::Environment walberlaEnv(argc, argv); + mpi::MPIManager::instance()->useWorldComm(); + + logging::configureLogging(walberlaEnv.config()); + + // read parameters + auto domainSetup = walberlaEnv.config()->getOneBlock("DomainSetup"); + auto parameters = walberlaEnv.config()->getOneBlock("Parameters"); + + const real_t omega = parameters.getParameter< real_t >("omega"); + const real_t inletVelocity = parameters.getParameter< real_t >("inletVelocity"); + const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10)) + uint_c(1); + const uint_t refinementDepth = parameters.getParameter< uint_t >("refinementDepth", uint_c(1)); + + auto loggingParameters = walberlaEnv.config()->getOneBlock("Logging"); + bool writeSetupForestAndReturn = loggingParameters.getParameter("writeSetupForestAndReturn", false); + + auto remainingTimeLoggerFrequency = + parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(3.0)); // in seconds + + auto flowSetup = std::make_shared< Channel >(refinementDepth); + + SetupBlockForest setupBfs; + WALBERLA_LOG_INFO_ON_ROOT("Generating SetupBlockForest...") + createSetupBlockForest(setupBfs, domainSetup, *flowSetup); + + // Create structured block forest + Vector3< uint_t > cellsPerBlock = domainSetup.getParameter< Vector3< uint_t > >("cellsPerBlock"); + + WALBERLA_LOG_INFO_ON_ROOT("Creating structured block forest...") + auto bfs = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()), setupBfs); + auto blocks = std::make_shared< StructuredBlockForest >(bfs, cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2]); + blocks->createCellBoundingBoxes(); + + if (writeSetupForestAndReturn) + { + WALBERLA_LOG_INFO_ON_ROOT("Writing SetupBlockForest to VTK file") + WALBERLA_ROOT_SECTION() { vtk::writeDomainDecomposition(blocks, "FreeSlipRefinementDomainDecomposition", "vtk_out"); } + + WALBERLA_LOG_INFO_ON_ROOT("Blocks created: " << setupBfs.getNumberOfBlocks()) + for (uint_t level = 0; level <= refinementDepth; level++) + { + WALBERLA_LOG_INFO_ON_ROOT("Level " << level << " Blocks: " << setupBfs.getNumberOfBlocks(level)) + } + WALBERLA_LOG_INFO_ON_ROOT("Ending program") + return EXIT_SUCCESS; + } + + WALBERLA_LOG_INFO_ON_ROOT("Blocks created: " << setupBfs.getNumberOfBlocks()) + for (uint_t level = 0; level <= refinementDepth; level++) + { + WALBERLA_LOG_INFO_ON_ROOT("Level " << level << " Blocks: " << setupBfs.getNumberOfBlocks(level)) + } + + StorageSpecification_T StorageSpec = StorageSpecification_T(); + BlockDataID pdfFieldId = lbm_generated::addPdfFieldToStorage(blocks, "pdf field", StorageSpec, uint_c(2)); + 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(3)); + + SweepCollection_T sweepCollection(blocks, pdfFieldId, velFieldId, omega); + for (auto& block : *blocks) + { + sweepCollection.initialise(&block); + } + + const FlagUID fluidFlagUID("Fluid"); + flowSetup->setupBoundaryFlagField(*blocks, flagFieldId); + geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID, 2); + BoundaryCollection_T boundaryCollection(blocks, flagFieldId, pdfFieldId, fluidFlagUID, inletVelocity); + + WALBERLA_LOG_INFO_ON_ROOT("Setting up communication...") + auto comm = + std::make_shared< blockforest::communication::NonUniformBufferedScheme< CommunicationStencil_T > >(blocks); + auto packInfo = lbm_generated::setupNonuniformPdfCommunication< PdfField_T >(blocks, pdfFieldId); + comm->addPackInfo(packInfo); + + lbm_generated::BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T > timestep( + blocks, pdfFieldId, sweepCollection, boundaryCollection, comm, packInfo); + + SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps); + uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0); + if (vtkWriteFrequency > 0) + { + auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "FreeSlipRefinementVTK", 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"); + } + timeloop.addFuncAfterTimeStep(timestep); + + // log remaining time + if (remainingTimeLoggerFrequency > 1.0) + { + timeloop.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency), + "remaining time logger"); + } + + WALBERLA_LOG_INFO_ON_ROOT("Starting Simulation with " << timesteps << " timesteps") + + timeloop.run(); + + + for (auto& block : *blocks) + { + sweepCollection.calculateMacroscopicParameters(&block); + Block& b = dynamic_cast< Block& >(block); + uint_t level = b.getLevel(); + + auto velField = b.getData< VectorField_T >(velFieldId); + for( auto it = velField->beginXYZ(); it != velField->end(); ++it ) + { + Cell localCell = it.cell(); + Cell globalCell(localCell); + blocks->transformBlockLocalToGlobalCell(globalCell, b); + + if (globalCell.y() >= (cell_idx_c(1 << level))) + { + WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(it.getF(0), inletVelocity, real_c(1e-5)); + } + } + } + return EXIT_SUCCESS; +} diff --git a/tests/lbm_generated/FreeSlipRefinement.prm b/tests/lbm_generated/FreeSlipRefinement.prm new file mode 100644 index 00000000..9149c34d --- /dev/null +++ b/tests/lbm_generated/FreeSlipRefinement.prm @@ -0,0 +1,26 @@ +Parameters +{ + omega 1.95; + inletVelocity 0.05; + timesteps 1000; + refinementDepth 1; + + remainingTimeLoggerFrequency 0; // in seconds + vtkWriteFrequency 0; +} + +DomainSetup +{ + domainSize <32, 16, 16>; + rootBlocks <4, 2, 2>; + + cellsPerBlock < 8, 8, 8 >; + periodic < 0, 0, 1 >; +} + +Logging +{ + logLevel info; // info progress detail tracing + writeSetupForestAndReturn false; +} + diff --git a/tests/lbm_generated/FreeSlipRefinement.py b/tests/lbm_generated/FreeSlipRefinement.py new file mode 100644 index 00000000..e695a8ae --- /dev/null +++ b/tests/lbm_generated/FreeSlipRefinement.py @@ -0,0 +1,49 @@ +import sympy as sp + +from pystencils import Target +from pystencils import fields + +from lbmpy.advanced_streaming.utility import get_timesteps +from lbmpy.boundaries import FreeSlip, UBB, ExtrapolationOutflow +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 + +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) + + freeslip = lbm_boundary_generator("FreeSlip", flag_uid="FreeSlip", boundary_object=FreeSlip(stencil)) + ubb = lbm_boundary_generator(class_name='UBB', flag_uid='UBB', + boundary_object=UBB([sp.Symbol("u_x"), 0.0, 0.0], data_type=data_type)) + outflow = lbm_boundary_generator(class_name='Outflow', flag_uid='Outflow', + boundary_object=ExtrapolationOutflow(stencil[4], method), + field_data_type=data_type) + + generate_lbm_package(ctx, name="FreeSlipRefinement", + collision_rule=collision_rule, + lbm_config=lbm_config, lbm_optimisation=lbm_opt, + nonuniform=True, boundaries=[freeslip, ubb, outflow], + macroscopic_fields=macroscopic_fields, + data_type=data_type, pdfs_data_type=data_type) + + generate_info_header(ctx, 'FreeSlipRefinementInfoHeader')