Forked from
waLBerla / waLBerla
983 commits behind the upstream repository.
-
Markus Holzer authoredMarkus Holzer authored
CodegenJacobiCPU.cpp 5.66 KiB
//======================================================================================================================
//
// 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 CodegenJacobiGPU.cpp
//! \author Martin Bauer <martin.bauer@fau.de>
//
//======================================================================================================================
#include "JacobiKernel2D.h"
#include "JacobiKernel3D.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 testJacobi2D()
{
uint_t xSize = 20;
uint_t ySize = 20;
// 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
true, true, true ); // full periodicity
BlockDataID fieldID = field::addToStorage<ScalarField>(blocks, "Field", real_t(0.0), field::fzyx);
// Initialize a quarter of the field with ones, the rest remains 0
// Jacobi averages the domain -> every cell should be at 0.25 at sufficiently many timesteps
for(auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
auto f = blockIt->getData<ScalarField>( fieldID );
for( cell_idx_t y = 0; y < cell_idx_c( f->ySize() / 2 ); ++y )
for( cell_idx_t x = 0; x < cell_idx_c( f->xSize() / 2 ); ++x )
f->get( x, y, 0 ) = 1.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(800);
SweepTimeloop timeloop ( blocks, numberOfTimesteps );
// Registering the sweep
timeloop.add() << BeforeFunction( commScheme, "Communication" )
<< 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));
}
void testJacobi3D()
{
uint_t xSize = 12;
uint_t ySize = 12;
uint_t zSize = 12;
// 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, zSize, // 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
true, true, true ); // no periodicity
BlockDataID fieldID = field::addToStorage<ScalarField>(blocks, "Field", real_t(0.0), field::fzyx);
// Initialize a quarter of the field with ones, the rest remains 0
// Jacobi averages the domain -> every cell should be at 0.25 at sufficiently many timesteps
for(auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
auto f = blockIt->getData<ScalarField>( fieldID );
for( cell_idx_t z = 0; z < cell_idx_c( f->zSize() / 2); ++z )
for( cell_idx_t y = 0; y < cell_idx_c( f->ySize() / 2 ); ++y )
for( cell_idx_t x = 0; x < cell_idx_c( f->xSize() / 2 ); ++x )
f->get( x, y, z ) = 1.0;
}
typedef blockforest::communication::UniformBufferedScheme<stencil::D3Q7> CommScheme;
typedef field::communication::PackInfo<ScalarField> Packing;
CommScheme commScheme(blocks);
commScheme.addDataToCommunicate( make_shared<Packing>(fieldID) );
// Create Timeloop
const uint_t numberOfTimesteps = uint_t(800); // number of timesteps for non-gui runs
SweepTimeloop timeloop ( blocks, numberOfTimesteps );
// Registering the sweep
timeloop.add() << BeforeFunction( commScheme, "Communication" )
<< Sweep( pystencils::JacobiKernel3D(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 / 8.0));
}
int main( int argc, char ** argv )
{
mpi::Environment env( argc, argv );
debug::enterTestMode();
testJacobi2D();
testJacobi3D();
return 0;
}