diff --git a/tests/field/CMakeLists.txt b/tests/field/CMakeLists.txt index bcb806a0d899c2e0dad75ed0f42b317bff5733e0..c0091537003048ed7f135181b7a96d55b7a3f15d 100644 --- a/tests/field/CMakeLists.txt +++ b/tests/field/CMakeLists.txt @@ -68,4 +68,9 @@ 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 CodegenPoissonGeneratedKernel FILE codegen/Poisson.py + OUT_FILES Poisson.cpp Poisson.h ) +waLBerla_compile_test( FILES codegen/CodegenPoisson.cpp DEPENDS gui timeloop CodegenPoissonGeneratedKernel) +waLBerla_execute_test( NAME CodegenPoisson ) + diff --git a/tests/field/codegen/CodegenJacobiCPU.cpp b/tests/field/codegen/CodegenJacobiCPU.cpp index aefaf67e2b0152621dad007e3a4c7dfa2422fcde..20bc8b06116b9437de17e5967bf44e83df4100fb 100644 --- a/tests/field/codegen/CodegenJacobiCPU.cpp +++ b/tests/field/codegen/CodegenJacobiCPU.cpp @@ -79,12 +79,13 @@ void testJacobi2D() // Registering the sweep timeloop.add() << BeforeFunction( commScheme, "Communication" ) - << Sweep( pystencils::JacobiKernel2D(fieldID, 1.0), "Jacobi Kernel" ); + << Sweep( pystencils::JacobiKernel2D(fieldID), "Jacobi Kernel" ); timeloop.run(); auto firstBlock = blocks->begin(); auto f = firstBlock->getData<ScalarField>( fieldID ); + WALBERLA_CHECK_FLOAT_EQUAL(f->get(0,0,0), real_t(1.0 / 4.0)); } @@ -127,7 +128,7 @@ void testJacobi3D() // Registering the sweep timeloop.add() << BeforeFunction( commScheme, "Communication" ) - << Sweep( pystencils::JacobiKernel3D(fieldID, 1.0), "Jacobi Kernel" ); + << Sweep( pystencils::JacobiKernel3D(fieldID), "Jacobi Kernel" ); timeloop.run(); diff --git a/tests/field/codegen/CodegenPoisson.cpp b/tests/field/codegen/CodegenPoisson.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c390831b649d68a075eb539333fdd8cfe3b4236f --- /dev/null +++ b/tests/field/codegen/CodegenPoisson.cpp @@ -0,0 +1,93 @@ +//====================================================================================================================== +// +// 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 CodegenPoisson.cpp +//! \author Markus Holzer <markus.holzer@fau.de> +// +//====================================================================================================================== + +#include "Poisson.h" +#include "blockforest/Initialization.h" +#include "blockforest/communication/UniformDirectScheme.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 "gui/Gui.h" + +#include "stencil/D2Q9.h" +#include "stencil/D3Q7.h" + +#include "timeloop/SweepTimeloop.h" + + +using namespace walberla; + +typedef GhostLayerField<double,1> ScalarField; + + +void testPoisson() +{ + uint_t L0 = 50; + uint_t xSize = L0; + uint_t ySize = L0; + + double h = 1 / (double(L0) + 1); + double rhs = 1.0; + // 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_t(1), // dx: length of one cell in physical coordinates + false, // one block per process - "false" means all blocks to one process + false, false, false ); // no periodicity + + + BlockDataID fieldID = field::addToStorage<ScalarField>(blocks, "Field", real_t(0.0)); + + typedef blockforest::communication::UniformBufferedScheme<stencil::D2Q9> CommScheme; + typedef field::communication::PackInfo<ScalarField> Packing; + CommScheme commScheme(blocks); + commScheme.addDataToCommunicate( make_shared<Packing>(fieldID) ); + + // Create Timeloop + const uint_t numberOfTimesteps = uint_t(2500); + SweepTimeloop timeloop ( blocks, numberOfTimesteps ); + + // Registering the sweep + timeloop.add() << BeforeFunction( commScheme, "Communication" ) + << Sweep( pystencils::Poisson(fieldID, h, rhs), "Poisson Kernel" ); + + timeloop.run(); + + auto firstBlock = blocks->begin(); + auto f = firstBlock->getData<ScalarField>( fieldID ); + WALBERLA_CHECK_FLOAT_EQUAL(f->get(int(L0/2), int(L0/2), 0), real_t(7.28886456002774547e-02)); +} + + +int main( int argc, char ** argv ) +{ + mpi::Environment env( argc, argv ); + debug::enterTestMode(); + + testPoisson(); + + return 0; +} diff --git a/tests/field/codegen/JacobiKernel.py b/tests/field/codegen/JacobiKernel.py index 4aa645e1d85e5664411a38561647b90ad03fc66e..c135ded28942cd897e39faed21b259b8073e2adb 100644 --- a/tests/field/codegen/JacobiKernel.py +++ b/tests/field/codegen/JacobiKernel.py @@ -1,26 +1,24 @@ -import sympy as sp +import numpy as np import pystencils as ps from pystencils_walberla import CodeGeneration, generate_sweep with CodeGeneration() as ctx: - h = sp.symbols("h") - - # ----- Jacobi 2D - created by specifying weights in nested list -------------------------- + # ----- Stencil 2D - created by specifying weights in nested list -------------------------- src, dst = ps.fields("src, src_tmp: [2D]", layout='fzyx') - stencil = [[0, 1, 0], - [1, 0, 1], - [0, 1, 0]] - assignments = ps.assignment_from_stencil(stencil, src, dst, normalization_factor=1 / (4 * h ** 2)) + stencil = [[1, 2, 3], + [4, 5, 6], + [7, 8, 9]] + assignments = ps.assignment_from_stencil(stencil, src, dst, normalization_factor=1 / np.sum(stencil)) generate_sweep(ctx, 'JacobiKernel2D', assignments, field_swaps=[(src, dst)]) - # ----- Jacobi 3D - created by using kernel_decorator with assignments in '@=' format ----- - src, dst = ps.fields("src, src_tmp: [3D]") + # ----- Stencil 3D - created by using kernel_decorator with assignments in '@=' format ----- + src, dst = ps.fields("src, src_tmp: [3D]", layout='fzyx') @ps.kernel def kernel_func(): - dst[0, 0, 0] @= (src[1, 0, 0] + src[-1, 0, 0] + - src[0, 1, 0] + src[0, -1, 0] + - src[0, 0, 1] + src[0, 0, -1]) / (6 * h ** 2) + dst[0, 0, 0] @= (3 * src[1, 0, 0] + 4 * src[-1, 0, 0] + + 5 * src[0, 1, 0] + 6 * src[0, -1, 0] + + 7 * src[0, 0, 1] + 8 * src[0, 0, -1]) / 33 generate_sweep(ctx, 'JacobiKernel3D', kernel_func, field_swaps=[(src, dst)]) diff --git a/tests/field/codegen/Poisson.py b/tests/field/codegen/Poisson.py new file mode 100644 index 0000000000000000000000000000000000000000..9cf604064eac2e8feff950e841642f9cf9b31e3d --- /dev/null +++ b/tests/field/codegen/Poisson.py @@ -0,0 +1,16 @@ +import sympy as sp +import pystencils as ps +from pystencils_walberla import CodeGeneration, generate_sweep + + +with CodeGeneration() as ctx: + # ----- Solving the 2D Poisson equation with constant rhs -------------------------- + rhs = sp.Symbol("rhs") + h = sp.Symbol("h") + src, dst = ps.fields("src, src_tmp: [2D]", layout='fzyx') + + @ps.kernel + def kernel_func(): + dst[0, 0] @= (src[1, 0] + src[-1, 0] + + src[0, 1] + src[0, -1]) / 4 + ((rhs * h**2) / 4) + generate_sweep(ctx, 'Poisson', kernel_func, field_swaps=[(src, dst)])