diff --git a/python/pystencils_walberla/__init__.py b/python/pystencils_walberla/__init__.py index f78f7fcf244e7fd140cd2abcc93ebaebaea2f94f..158a294523f207daeed3b977408a08a2157ece72 100644 --- a/python/pystencils_walberla/__init__.py +++ b/python/pystencils_walberla/__init__.py @@ -1,4 +1,4 @@ -from .boundary import generate_staggered_boundary, generate_staggered_flux_boundary +from .boundary import generate_boundary, generate_staggered_boundary, generate_staggered_flux_boundary from .cmake_integration import CodeGeneration, ManualCodeGenerationContext from .function_generator import function_generator @@ -8,7 +8,7 @@ from .pack_info import (generate_pack_info, generate_pack_info_for_field, generate_pack_info_from_kernel, generate_mpidtype_info_from_kernel) from .utility import generate_info_header, get_vectorize_instruction_set, config_from_context -__all__ = ['generate_staggered_boundary', 'generate_staggered_flux_boundary', +__all__ = ['generate_boundary', 'generate_staggered_boundary', 'generate_staggered_flux_boundary', 'CodeGeneration', 'ManualCodeGenerationContext', 'function_generator', 'generate_sweep', 'generate_selective_sweep', 'generate_sweep_collection', diff --git a/tests/field/CMakeLists.txt b/tests/field/CMakeLists.txt index b48f4ac79d1a778ccd1cadad910b6fab00a99b36..f2402a3e4e3a4aedc33184a43fcfba21d36cb6eb 100644 --- a/tests/field/CMakeLists.txt +++ b/tests/field/CMakeLists.txt @@ -71,6 +71,11 @@ waLBerla_generate_target_from_python(NAME CodegenJacobiCPUGeneratedJacobiKernel waLBerla_compile_test( FILES codegen/CodegenJacobiCPU.cpp DEPENDS gui timeloop CodegenJacobiCPUGeneratedJacobiKernel) waLBerla_execute_test( NAME CodegenJacobiCPU ) +waLBerla_generate_target_from_python(NAME CodegenADEGenerated FILE codegen/ADEKernel.py + OUT_FILES ADE.cpp ADE.h Dirichlet.cpp Dirichlet.h Neumann.cpp Neumann.h ADEPackInfo.cpp ADEPackInfo.h) +waLBerla_compile_test( FILES codegen/ADE.cpp DEPENDS blockforest core field stencil timeloop vtk CodegenADEGenerated) +waLBerla_execute_test( NAME CodegenADE ) + waLBerla_generate_target_from_python(NAME SweepCollectionKernel FILE codegen/SweepCollection.py OUT_FILES SweepCollection.h SweepCollection.cpp) waLBerla_compile_test( FILES codegen/SweepCollection.cpp DEPENDS timeloop SweepCollectionKernel) diff --git a/tests/field/codegen/ADE.cpp b/tests/field/codegen/ADE.cpp new file mode 100644 index 0000000000000000000000000000000000000000..29bada362646851e2f1230f6501b2746954880f5 --- /dev/null +++ b/tests/field/codegen/ADE.cpp @@ -0,0 +1,187 @@ +//====================================================================================================================== +// +// 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 ADE.cpp +//! \author Markus Holzer <markus.holzer@fau.de> +// +//====================================================================================================================== +#include "blockforest/Initialization.h" +#include "blockforest/communication/UniformBufferedScheme.h" + +#include "core/Environment.h" +#include "core/debug/TestSubsystem.h" + +#include "field/AddToStorage.h" +#include "field/communication/PackInfo.h" +#include "field/vtk/VTKWriter.h" + +#include "stencil/D2Q5.h" + +#include "timeloop/SweepTimeloop.h" + +#include "vtk/VTKOutput.h" + +#include "ADE.h" +#include "Dirichlet.h" +#include "Neumann.h" +#include "ADEPackInfo.h" + + +using namespace walberla; + +using ScalarField = GhostLayerField<real_t, 1>; +using VectorField = GhostLayerField<real_t, 2>; +using flag_t = walberla::uint8_t; +using FlagField_T = FlagField< flag_t >; + +using CommScheme = blockforest::communication::UniformBufferedScheme<stencil::D2Q5>; + +void setupBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID, + const FlagUID& fluidFlagUID, const FlagUID& dirichletZeroFlagUID, const FlagUID& dirichletOneFlagUID, const FlagUID& neumannFlagUID) +{ + + const AABB & domain = sbfs->getDomain(); + const AABB inflowAABB( 10, 25, 0, 30, 39, domain.zMax()); + + for (auto bIt = sbfs->begin(); bIt != sbfs->end(); ++bIt) + { + + Block& b = dynamic_cast< Block& >(*bIt); + auto flagField = b.getData< FlagField_T >(flagFieldID); + + uint8_t fluidFlag = flagField->registerFlag(fluidFlagUID); + uint8_t dirichletZeroFlag = flagField->registerFlag(dirichletZeroFlagUID); + uint8_t dirichletOneFlag = flagField->registerFlag(dirichletOneFlagUID); + uint8_t neumannFlag = flagField->registerFlag(neumannFlagUID); + + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, { + Cell globalCell = Cell(x, y, z); + sbfs->transformBlockLocalToGlobalCell(globalCell, *bIt); + const Vector3< real_t > globalPoint = sbfs->getCellCenter(globalCell); + + if(globalPoint[0] < 0 || globalPoint[1] < 0 || globalPoint[1] > domain.yMax()) + { + addFlag(flagField->get(x, y, z), dirichletZeroFlag); + } + else if (globalPoint[0] > domain.xMax()) + { + addFlag(flagField->get(x, y, z), neumannFlag); + } + else if (inflowAABB.contains(globalPoint)) + { + addFlag(flagField->get(x, y, z), dirichletOneFlag); + } + else + { + addFlag(flagField->get(x, y, z), fluidFlag); + } + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ + + } +} + +int main( int argc, char ** argv ) +{ + mpi::Environment env( argc, argv ); + debug::enterTestMode(); + + logging::Logging::instance()->setLogLevel( logging::Logging::INFO ); + + uint_t xSize = 128; + uint_t ySize = 64; + // Create blocks + shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid ( + uint_t(1) , uint_t(1), uint_t(1), // number of blocks in x,y,z direction + xSize, ySize, uint_t(1), // how many cells per block (x,y,z) + real_c(1.0), // dx: length of one cell in physical coordinates + false, // one block per process - "false" means all blocks to one process + false, false, false ); // full periodicity + + real_t diffusivity = real_c(0.01); + real_t dx = real_c(1.0); + real_t dt = real_c(0.01); + + BlockDataID src = field::addToStorage<ScalarField>(blocks, "src", real_c(0.0)); + BlockDataID vel = field::addToStorage<VectorField>(blocks, "vel", real_c(0.0)); + BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field"); + + for (auto& block : *blocks) + { + const AABB & domain = blocks->getDomain(); + const AABB inflowAABB( 10, 25, 0, 30, 39, domain.zMax()); + + auto * velField = block.getData< VectorField > ( vel ); + auto * srcField = block.getData< ScalarField > ( src ); + + for( auto it = velField->beginWithGhostLayer(1); it != velField->end(); ++it ) + { + velField->get(it.cell(), 0) = real_c(1.0); + if (inflowAABB.contains( it.x(), it.y(), it.z() )) + { + srcField->get(it.cell()) = real_c(1.0); + } + } + } + + const FlagUID fluidFlagUID("Fluid"); + const FlagUID dirichletZeroFlagUID("dirichletZero"); + const FlagUID dirichletOneFlagUID("dirichletOne"); + const FlagUID neumannFlagUID("neumann"); + + setupBoundary(blocks, flagFieldID, fluidFlagUID, dirichletZeroFlagUID, dirichletOneFlagUID, neumannFlagUID); + + pystencils::Dirichlet dirichletZero(blocks, src, real_c(0.0)); + pystencils::Dirichlet dirichletOne(blocks, src, real_c(1.0)); + pystencils::Neumann neumann(blocks, src); + + dirichletZero.fillFromFlagField< FlagField_T >(blocks, flagFieldID, dirichletZeroFlagUID, fluidFlagUID); + dirichletOne.fillFromFlagField< FlagField_T >(blocks, flagFieldID, dirichletOneFlagUID, fluidFlagUID); + neumann.fillFromFlagField< FlagField_T >(blocks, flagFieldID, neumannFlagUID, fluidFlagUID); + + CommScheme commScheme(blocks); + commScheme.addDataToCommunicate( make_shared<pystencils::ADEPackInfo>(src, vel) ); + + // Create Timeloop + const uint_t numberOfTimesteps = uint_t(100); // number of timesteps for non-gui runs + SweepTimeloop timeloop ( blocks, numberOfTimesteps ); + + // Registering the sweep + timeloop.add() << BeforeFunction( commScheme, "Communication" ) + << Sweep(dirichletZero, "dirichlet zero"); + timeloop.add() << Sweep(dirichletOne, "dirichlet one" ); + timeloop.add() << Sweep(neumann, "neumann" ); + timeloop.add() << Sweep( pystencils::ADE(src, vel, diffusivity, dt, dx), "advection diffusion kernel" ); + + const uint_t vtkWriteFrequency = 0; + if (vtkWriteFrequency > 0) + { + auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false, + "vtk_ADE", "simulation_step", false, true, true, + false, 0); + + auto concentrationWriter = make_shared< field::VTKWriter< ScalarField, float32 > >(src, "concentration"); + vtkOutput->addCellDataWriter(concentrationWriter); + + auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(flagFieldID, "flag"); + vtkOutput->addCellDataWriter(flagWriter); + + + timeloop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output"); + } + timeloop.run(); + + + return EXIT_SUCCESS; +} diff --git a/tests/field/codegen/ADEKernel.py b/tests/field/codegen/ADEKernel.py new file mode 100644 index 0000000000000000000000000000000000000000..60da9cd63594084b897e9bd5f0443cffc7bf621e --- /dev/null +++ b/tests/field/codegen/ADEKernel.py @@ -0,0 +1,25 @@ +import sympy as sp +import pystencils as ps +import pystencils_walberla as psw + +with psw.CodeGeneration() as ctx: + dim = 2 + c, c_tmp, u = ps.fields(f"c, c_tmp, u({dim}): [{dim}d]") + D, dx, dt = sp.symbols("D dx dt") + + pde = ps.fd.transient(c) - ps.fd.diffusion(c, D) + ps.fd.advection(c, u) + + discretize = ps.fd.Discretization2ndOrder(dx, dt) + assign = [ps.Assignment(c_tmp.center(), discretize(pde)),] + + target = ps.Target.GPU if ctx.gpu else ps.Target.CPU + + psw.generate_sweep(ctx, 'ADE', assign, field_swaps=[(c, c_tmp)], target=target) + + dirichlet = ps.Dirichlet(sp.Symbol("concentration")) + neumann = ps.Neumann() + + psw.generate_boundary(ctx, 'Dirichlet', dirichlet, field=c, target=target) + psw.generate_boundary(ctx, 'Neumann', neumann, field=c, target=target) + + psw.generate_pack_info_from_kernel(ctx, 'ADEPackInfo', assign, target=target) \ No newline at end of file