added Poisson solver testcase

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 )
+
 

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();
 

tests/field/codegen/CodegenPoisson.cpp

+//======================================================================================================================
+//
+//  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;
+}

tests/field/codegen/JacobiKernel.py

-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)])

tests/field/codegen/Poisson.py

+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)])