diff --git a/src/pe/rigidbody/Squirmer.cpp b/src/pe/rigidbody/Squirmer.cpp
index 9bd8098023243c39f1ca40c9afb9d85be95b821d..99f650b00290b619fa411c1a74f77a8c193c6200 100644
--- a/src/pe/rigidbody/Squirmer.cpp
+++ b/src/pe/rigidbody/Squirmer.cpp
@@ -106,7 +106,7 @@ const Vec3 Squirmer::velFromBF( real_t px, real_t py, real_t pz ) const
*/
const Vec3 Squirmer::velFromBF( const Vec3& rpos ) const
{
- return Sphere::velFromBF( rpos ) + getSquirmerVelocity( rpos );
+ return Sphere::velFromBF( rpos ) + getSquirmerVelocity( R_ * rpos );
}
//*************************************************************************************************
@@ -130,6 +130,7 @@ const Vec3 Squirmer::velFromWF( real_t px, real_t py, real_t pz ) const
//*************************************************************************************************
/*!\brief Returns the squirmer's surface velocity at a given relative position.
*
+ * \param rpos The relative global coordinate.
* \return The local surface velocity of the squirmer.
*/
const Vec3 Squirmer::getSquirmerVelocity( const Vec3& rpos ) const
@@ -153,8 +154,7 @@ const Vec3 Squirmer::getSquirmerVelocity( const Vec3& rpos ) const
*/
const Vec3 Squirmer::velFromWF( const Vec3& gpos ) const
{
- const auto rpos = pointFromWFtoBF( gpos );
- return Sphere::velFromWF( gpos ) + getSquirmerVelocity( rpos );
+ return Sphere::velFromWF( gpos ) + getSquirmerVelocity( gpos - gpos_ );
}
//*************************************************************************************************
diff --git a/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h b/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h
index a6807a6ff0a78df2cf23f7831a298821cd1bb0eb..161dfc872e1094df917d2394048790afdf52da92 100644
--- a/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h
+++ b/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h
@@ -59,8 +59,6 @@ namespace pe_coupling {
template< typename LatticeModel_T, typename FlagField_T >
class SimpleBB : public Boundary< typename FlagField_T::flag_t >
{
- static_assert( LatticeModel_T::compressible == false, "Only works with incompressible models!" );
-
typedef lbm::PdfField< LatticeModel_T > PDFField_T;
typedef typename LatticeModel_T::Stencil Stencil_T;
typedef typename FlagField_T::flag_t flag_t;
@@ -172,9 +170,21 @@ inline void SimpleBB< LatticeModel_T, FlagField_T >::treatDirection( const cell_
auto boundaryVelocity = body.velFromWF( boundaryPosition );
// include effect of boundary velocity
- pdf_new -= real_c(3.0) * alpha * LatticeModel_T::w[ Stencil_T::idx[dir] ] * ( real_c( stencil::cx[ dir ] ) * boundaryVelocity[0] +
- real_c( stencil::cy[ dir ] ) * boundaryVelocity[1] +
- real_c( stencil::cz[ dir ] ) * boundaryVelocity[2] );
+ if( LatticeModel_T::compressible )
+ {
+ const auto density = pdfField_->getDensity(x,y,z);
+ pdf_new -= real_c(3.0) * alpha * density * LatticeModel_T::w[ Stencil_T::idx[dir] ] *
+ ( real_c( stencil::cx[ dir ] ) * boundaryVelocity[0] +
+ real_c( stencil::cy[ dir ] ) * boundaryVelocity[1] +
+ real_c( stencil::cz[ dir ] ) * boundaryVelocity[2] );
+ }
+ else
+ {
+ pdf_new -= real_c(3.0) * alpha * LatticeModel_T::w[ Stencil_T::idx[dir] ] *
+ ( real_c( stencil::cx[ dir ] ) * boundaryVelocity[0] +
+ real_c( stencil::cy[ dir ] ) * boundaryVelocity[1] +
+ real_c( stencil::cz[ dir ] ) * boundaryVelocity[2] );
+ }
// carry out the boundary handling
pdfField_->get( nx, ny, nz, Stencil_T::invDirIdx(dir) ) = pdf_new;
diff --git a/tests/pe_coupling/CMakeLists.txt b/tests/pe_coupling/CMakeLists.txt
index 24b8d2a53791b812d2aa51b363e60408357103dc..e31ded9e826492cb2926ada9d7272dd4647bccd6 100644
--- a/tests/pe_coupling/CMakeLists.txt
+++ b/tests/pe_coupling/CMakeLists.txt
@@ -123,3 +123,6 @@ waLBerla_execute_test( NAME TorqueSphereMEMMRFuncTest COMMAND $<TARGET_FILE
waLBerla_execute_test( NAME TorqueSphereMEMMRSingleTest COMMAND $<TARGET_FILE:TorqueSphereMEM> --MO_MR PROCESSES 1 LABELS longrun CONFIGURATIONS Release RelWithDbgInfo )
waLBerla_execute_test( NAME TorqueSphereMEMMRParallelTest COMMAND $<TARGET_FILE:TorqueSphereMEM> --MO_MR PROCESSES 8 LABELS verylongrun CONFIGURATIONS Release RelWithDbgInfo )
+waLBerla_compile_test( FILES momentum_exchange_method/SquirmerTest.cpp DEPENDS blockforest pe timeloop )
+waLBerla_execute_test( NAME SquirmerShortTest COMMAND $<TARGET_FILE:SquirmerTest> --shortrun PROCESSES 1 )
+waLBerla_execute_test( NAME SquirmerTest COMMAND $<TARGET_FILE:SquirmerTest> PROCESSES 1 CONFIGURATIONS Release RelWithDbgInfo )
diff --git a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..235fcc2a62739f4d5cb559e164a34cf143db82a3
--- /dev/null
+++ b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
@@ -0,0 +1,474 @@
+//======================================================================================================================
+//
+// 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 SquirmerTest.cpp
+//! \ingroup pe_coupling
+//! \author Christian Burkard <buch@icp.uni-stuttgart.de>
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/all.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/Debug.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/logging/Logging.h"
+#include "core/math/all.h"
+#include "core/SharedFunctor.h"
+#include "core/timing/RemainingTimeLogger.h"
+#include "core/mpi/MPIManager.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/AddToStorage.h"
+#include "field/communication/PackInfo.h"
+
+#include "lbm/boundary/all.h"
+#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/lattice_model/ForceModel.h"
+#include "lbm/sweeps/CellwiseSweep.h"
+#include "lbm/sweeps/SweepWrappers.h"
+
+#include "stencil/D3Q7.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "pe/rigidbody/SquirmerFactory.h"
+#include "pe/communication/rigidbody/Squirmer.h"
+#include "pe/basic.h"
+
+#include "pe_coupling/mapping/all.h"
+#include "pe_coupling/momentum_exchange_method/all.h"
+#include "pe_coupling/utility/all.h"
+
+#include "field/vtk/FlagFieldCellFilter.h"
+#include "field/vtk/VTKWriter.h"
+#include "lbm/vtk/Velocity.h"
+#include "vtk/VTKOutput.h"
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+using walberla::uint_t;
+
+// PDF field, flag field & body field
+typedef lbm::force_model::None ForceModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT, false, ForceModel_T > LatticeModel_T;
+
+typedef LatticeModel_T::Stencil Stencil_T;
+typedef lbm::PdfField< LatticeModel_T > PdfField_T;
+
+typedef walberla::uint8_t flag_t;
+typedef FlagField< flag_t > FlagField_T;
+typedef GhostLayerField< pe::BodyID, 1 > PeBodyField_T;
+
+const uint_t FieldGhostLayers = 1;
+
+// boundary handling
+typedef pe_coupling::SimpleBB < LatticeModel_T, FlagField_T > MO_BB_T;
+
+typedef boost::tuples::tuple< MO_BB_T > BoundaryConditions_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+
+typedef boost::tuple<pe::Squirmer> BodyTypeTuple ;
+
+///////////
+// FLAGS //
+///////////
+
+const FlagUID Fluid_Flag ( "fluid" );
+const FlagUID MO_BB_Flag ( "moving obstacle BB" );
+const FlagUID FormerMO_BB_Flag ( "former moving obstacle BB" );
+
+/////////////////////////////////////
+// BOUNDARY HANDLING CUSTOMIZATION //
+/////////////////////////////////////
+
+class MyBoundaryHandling
+{
+public:
+
+ MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID, const BlockDataID & bodyFieldID ) :
+ flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ), bodyFieldID_ ( bodyFieldID ) {}
+
+ BoundaryHandling_T * operator()( IBlock* const block, const StructuredBlockStorage* const storage ) const;
+
+private:
+
+ const BlockDataID flagFieldID_;
+ const BlockDataID pdfFieldID_;
+ const BlockDataID bodyFieldID_;
+
+}; // class MyBoundaryHandling
+
+////////////////////////////////
+// VTK OUTPUT HELPER FUNCTION //
+////////////////////////////////
+
+BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const StructuredBlockStorage * const storage ) const
+{
+ WALBERLA_ASSERT_NOT_NULLPTR( block );
+ WALBERLA_ASSERT_NOT_NULLPTR( storage );
+
+ FlagField_T * flagField = block->getData< FlagField_T >( flagFieldID_ );
+ PdfField_T * pdfField = block->getData< PdfField_T > ( pdfFieldID_ );
+ PeBodyField_T * bodyField = block->getData< PeBodyField_T >( bodyFieldID_ );
+
+ const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
+
+ BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
+ boost::tuples::make_tuple( MO_BB_T ( "MO_BB", MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+
+ handling->fillWithDomain( FieldGhostLayers );
+
+ return handling;
+}
+
+shared_ptr< vtk::VTKOutput> createFluidFieldVTKWriter( shared_ptr< StructuredBlockForest > & blocks,
+ const BlockDataID & pdfFieldId, const BlockDataID & flagFieldId,
+ std::string fname )
+{
+ auto pdfFieldVTKWriter = vtk::createVTKOutput_BlockData( blocks, fname, uint_c(1), uint_c(0) );
+
+ field::FlagFieldCellFilter< FlagField_T > fluidFilter( flagFieldId );
+ fluidFilter.addFlag( Fluid_Flag );
+ pdfFieldVTKWriter->addCellInclusionFilter( fluidFilter );
+
+ auto velocityWriter = make_shared< lbm::VelocityVTKWriter< LatticeModel_T, float > >( pdfFieldId, "VelocityFromPDF" );
+ pdfFieldVTKWriter->addCellDataWriter( velocityWriter );
+
+ return pdfFieldVTKWriter;
+}
+
+class ResetPosition
+{
+public:
+
+ explicit ResetPosition( const shared_ptr<StructuredBlockStorage> & blockStorage,
+ const BlockDataID & bodyStorageID,
+ real_t dt )
+ : blockStorage_( blockStorage )
+ , bodyStorageID_( bodyStorageID )
+ , dt_( dt )
+ {}
+
+ void operator()()
+ {
+ for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
+ {
+ for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->setPosition( bodyIt->getPosition() - bodyIt->getLinearVel() * dt_ );
+ }
+ }
+ }
+
+private:
+
+ shared_ptr<StructuredBlockStorage> blockStorage_;
+ const BlockDataID & bodyStorageID_;
+
+ real_t dt_;
+
+}; // class ResetPosition
+
+
+Vector3<real_t> analytical_velocity(Vector3<real_t> r, real_t b1, real_t b2, Vector3<real_t> e, real_t radius)
+{
+ auto r_len = r.length();
+ auto rs = r.getNormalized();
+ auto frac = radius/r_len;
+ auto frac2 = frac*frac;
+ auto frac3 = frac2*frac;
+ auto frac4 = frac3*frac;
+
+ return b1 * frac3 * ((e*rs) * rs - e / real_c(3.0)) + ( frac4 - frac2 ) * b2 * 0.5 * ( 3 * (e * rs) * (e * rs) - 1 ) * rs + frac4 * b2 * (e*rs) * ( (e*rs) * rs - e);
+}
+
+
+int main( int argc, char **argv )
+{
+ debug::enterTestMode();
+
+ mpi::Environment env( argc, argv );
+
+ auto processes = MPIManager::instance()->numProcesses();
+
+ if( processes != 1 && processes != 2 && processes != 4 && processes != 8)
+ {
+ std::cerr << "Number of processes must be equal to either 1, 2, 4, or 8!" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ bool shortrun = false;
+ bool writevtk = false;
+ for( int i = 1; i < argc; ++i )
+ {
+ if( std::strcmp( argv[i], "--shortrun" ) == 0 ) { shortrun = true; continue; }
+ if( std::strcmp( argv[i], "--vtk" ) == 0 ) { writevtk = true; continue; }
+ }
+
+ ///////////////////////////
+ // SIMULATION PROPERTIES //
+ ///////////////////////////
+
+ real_t R = 3;
+ uint_t L = 24;
+ uint_t T = 1000;
+ if (shortrun)
+ {
+ L = 10;
+ T = 10;
+ }
+ const real_t dx = real_c(1.0);
+ const real_t dt = real_c(1.0);
+ const real_t omega = lbm::collision_model::omegaFromViscosity( real_c(1.0/6.0) );
+ const real_t v0 = real_c(0.01);
+ const real_t beta = real_c(5.0);
+
+ ///////////////////////////
+ // BLOCK STRUCTURE SETUP //
+ ///////////////////////////
+ const uint_t XBlocks = (processes >= 2) ? uint_t( 2 ) : uint_t( 1 );
+ const uint_t YBlocks = (processes >= 4) ? uint_t( 2 ) : uint_t( 1 );
+ const uint_t ZBlocks = (processes == 8) ? uint_t( 2 ) : uint_t( 1 );
+ const uint_t XCells = L / XBlocks;
+ const uint_t YCells = L / YBlocks;
+ const uint_t ZCells = L / ZBlocks;
+
+ // create fully periodic domain
+ auto blocks = blockforest::createUniformBlockGrid( XBlocks, YBlocks, ZBlocks, XCells, YCells, ZCells, dx, true,
+ true, true, true );
+
+ ////////
+ // PE //
+ ////////
+
+ auto globalBodyStorage = make_shared<pe::BodyStorage>();
+ pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
+ auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
+ auto ccdID = blocks->addBlockData(pe::ccd::createHashGridsDataHandling( globalBodyStorage, bodyStorageID ), "CCD");
+ auto fcdID = blocks->addBlockData(pe::fcd::createSimpleFCDDataHandling<BodyTypeTuple>(), "FCD");
+
+ pe::cr::HCSITS cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID);
+
+ /////////////////
+ // PE COUPLING //
+ /////////////////
+
+ // connect to pe
+ const real_t overlap = real_c( 1.5 ) * dx;
+
+ if( R > real_c( L ) * real_c(0.5) - overlap )
+ {
+ std::cerr << "Periodic sphere is too large and would lead to invalid PE state!" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ boost::function<void(void)> syncCall = boost::bind( pe::syncShadowOwners<BodyTypeTuple>, boost::ref(blocks->getBlockForest()), bodyStorageID, static_cast<WcTimingTree*>(NULL), overlap, false );
+
+ const auto myMat = pe::createMaterial( "myMat", real_c(1), real_t(0), real_t(1), real_t(1), real_t(0), real_t(1), real_t(1), real_t(0), real_t(0) );
+
+ // create the squirmer in the middle of the domain
+ const Vector3<real_t> position (real_c(L)*real_c(0.5), real_c(L)*real_c(0.5), real_c(L)*real_c(0.5));
+ const Vector3<real_t> up(0.0, 0.0, 1.0);
+ const Vector3<real_t> orientation(1.0, 0.0, 0.0);
+ const auto w = std::acos(up*orientation);
+ math::Quaternion<real_t> q(up % orientation, w);
+ auto squirmer = pe::createSquirmer(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, position, R, v0, beta, myMat);
+ if (squirmer)
+ squirmer->setOrientation(q);
+
+ syncCall();
+
+ ///////////////////////
+ // ADD DATA TO BLOCKS //
+ ////////////////////////
+
+ // create the lattice model
+ LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ) );
+
+ // add PDF field ( uInit = <0,0,0>, rhoInit = 1 )
+ BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel, FieldGhostLayers );
+
+ // add flag field
+ BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
+
+ // add body field
+ BlockDataID bodyFieldID = field::addToStorage< PeBodyField_T >( blocks, "body field", NULL, field::zyxf );
+
+ // add boundary handling & initialize outer domain boundaries
+ BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >(
+ MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
+
+ ///////////////
+ // TIME LOOP //
+ ///////////////
+
+ // create the timeloop
+ SweepTimeloop timeloop( blocks->getBlockStorage(), T );
+
+ // sweep for updating the pe body mapping into the LBM simulation
+ timeloop.add()
+ << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag, FormerMO_BB_Flag ),
+ "Body Mapping" );
+
+ // sweep for restoring PDFs in cells previously occupied by pe bodies
+ typedef pe_coupling::EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T > Reconstructor_T;
+ Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID );
+ timeloop.add()
+ << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,
+ reconstructor, FormerMO_BB_Flag, Fluid_Flag ), "PDF Restore" );
+
+ // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
+ boost::function< void () > commFunction;
+
+ blockforest::communication::UniformBufferedScheme< stencil::D3Q27 > scheme( blocks );
+ scheme.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldID ) );
+ commFunction = scheme;
+
+ // uses standard bounce back boundary conditions
+ pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
+
+ auto sweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+
+ // collision sweep
+ timeloop.add() << Sweep( lbm::makeCollideSweep( sweep ), "cell-wise LB sweep (collide)" );
+
+ // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
+ timeloop.add() << BeforeFunction( commFunction, "LBM Communication" )
+ << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
+
+ // streaming
+ timeloop.add() << Sweep( lbm::makeStreamSweep( sweep ), "cell-wise LB sweep (stream)" );
+
+ // add pe timesteps
+ timeloop.addFuncAfterTimeStep( pe_coupling::TimeStep( blocks, bodyStorageID, cr, syncCall, dt ), "pe Time Step" );
+ timeloop.addFuncAfterTimeStep( ResetPosition( blocks, bodyStorageID, dt ), "Reset pe positions" );
+
+ timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
+
+ ////////////////////////
+ // EXECUTE SIMULATION //
+ ////////////////////////
+
+ WALBERLA_LOG_INFO_ON_ROOT("Starting test...");
+ WcTimingPool timeloopTiming;
+ timeloop.run( timeloopTiming );
+ timeloopTiming.logResultOnRoot();
+
+ if (writevtk)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Writing out VTK file...");
+ auto pdfFieldVTKWriter = createFluidFieldVTKWriter( blocks, pdfFieldID, flagFieldID, "fluid_field" );
+ auto wf = vtk::writeFiles( pdfFieldVTKWriter );
+ wf();
+ auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", 1, 0 );
+ flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
+ vtk::writeFiles( flagFieldVTK )();
+ }
+
+ if (squirmer) // this object only exists on the node that contains position
+ {
+ WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getSquirmerVelocity(), v0);
+ WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getSquirmerBeta(), beta);
+ WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getRadius(), R);
+ WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getPosition(), position);
+ WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getQuaternion(), q);
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("Checking result...");
+ auto b1 = real_c(1.5) * v0;
+ auto b2 = beta * b1;
+ auto e = q.rotate(up).getNormalized();
+ auto radius = R;
+ auto squirmer_pos = position;
+
+ real_t abs_tolerance = real_c(0.0026);
+ real_t rel_tolerance = real_c(0.10);
+
+ for ( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
+ {
+ IBlock & block = *blockIt;
+ auto src = block.getData< PdfField_T >(pdfFieldID);
+
+ auto flagField = block.getData<FlagField_T> (flagFieldID);
+ auto flag = flagField->getFlag(Fluid_Flag);
+
+ WALBERLA_FOR_ALL_CELLS_XYZ(src, {
+ Vector3<real_t> pos;
+ blocks->getBlockLocalCellCenter( block, Cell(x,y,z), pos[0], pos[1], pos[2] );
+
+ if( flagField->isFlagSet(x, y, z, flag) ) //check for fluid/non-boundary
+ {
+ auto r = pos - squirmer_pos;
+
+ auto v_ana = analytical_velocity(r, b1, b2, e, radius);
+
+ // Superposition with the nearest neighbours (faces)
+ for( auto s = stencil::D3Q7::beginNoCenter(); s != stencil::D3Q7::end(); ++s )
+ {
+ Vector3<real_t> imidpoint(squirmer_pos);
+ imidpoint[0] += real_c(s.cx() * int_c(L));
+ imidpoint[1] += real_c(s.cy() * int_c(L));
+ imidpoint[2] += real_c(s.cz() * int_c(L));
+ v_ana += analytical_velocity( pos - imidpoint, b1, b2, e, radius );
+ }
+
+ if( !shortrun && r.length() > 2*radius ) // Exclude the volume around the squirmer, as the analytical solution is a far field solution only
+ {
+ auto v_data = src->getVelocity(x, y, z);
+ auto diff = v_data - v_ana;
+
+ for( uint_t d = 0; d < 3; ++d )
+ {
+ if( std::abs( diff[d] / v_ana[d] ) > rel_tolerance && std::abs( diff[d] ) > abs_tolerance )
+ {
+ WALBERLA_LOG_DEVEL("Difference too large in " << Cell(x,y,z) << " (" << r.length() << "). Expected " << v_ana << ", got " << v_data << ", relative " << std::abs(diff[d] / v_ana[d]) << ", absolute " << std::abs(diff[d]));
+ }
+ }
+ }
+
+ if( writevtk )
+ {
+ // store the value into the PDF field so we can write it to VTK later
+ src->setToEquilibrium( x, y, z, v_ana );
+ }
+ }
+ });
+ }
+
+ if (writevtk)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Writing out reference VTK file...");
+ auto pdfFieldVTKWriter = createFluidFieldVTKWriter( blocks, pdfFieldID, flagFieldID, "fluid_field_ref" );
+ auto wf = vtk::writeFiles( pdfFieldVTKWriter );
+ wf();
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("Completed test.");
+
+ return 0;
+}