apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp

@@ -16,6 +16,7 @@
 //! \file ComplexGeometry.cpp
 //! \ingroup mesh
 //! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//! \author Markus Holzer <markus.holzer@fau.de>
 //
 //======================================================================================================================
 
@@ -84,181 +85,208 @@
 #include "mesh_common/vtk/CommonDataSources.h"
 #include "mesh_common/vtk/VTKMeshWriter.h"
 
-namespace walberla {
+namespace walberla
+{
 
 template< typename MeshType >
-void vertexToFaceColor( MeshType & mesh, const typename MeshType::Color & defaultColor )
+void vertexToFaceColor(MeshType& mesh, const typename MeshType::Color& defaultColor)
 {
-   WALBERLA_CHECK( mesh.has_vertex_colors() );
+   WALBERLA_CHECK(mesh.has_vertex_colors())
    mesh.request_face_colors();
 
-   for( auto faceIt = mesh.faces_begin(); faceIt != mesh.faces_end(); ++faceIt )
+   for (auto faceIt = mesh.faces_begin(); faceIt != mesh.faces_end(); ++faceIt)
    {
       typename MeshType::Color vertexColor;
 
       bool useVertexColor = true;
 
-      auto vertexIt = mesh.fv_iter( *faceIt );
-      WALBERLA_ASSERT( vertexIt.is_valid() );
-      
-      vertexColor = mesh.color( *vertexIt );
+      auto vertexIt = mesh.fv_iter(*faceIt);
+      WALBERLA_ASSERT(vertexIt.is_valid())
+
+      vertexColor = mesh.color(*vertexIt);
 
       ++vertexIt;
-      while( vertexIt.is_valid() && useVertexColor )
+      while (vertexIt.is_valid() && useVertexColor)
       {
-         if( vertexColor != mesh.color( *vertexIt ) )
-            useVertexColor = false;
+         if (vertexColor != mesh.color(*vertexIt)) useVertexColor = false;
          ++vertexIt;
       }
 
-      mesh.set_color( *faceIt, useVertexColor ? vertexColor : defaultColor );
+      mesh.set_color(*faceIt, useVertexColor ? vertexColor : defaultColor);
    }
 }
 
-
-int main( int argc, char * argv[] )
+int main(int argc, char* argv[])
 {
-   Environment env( argc, argv );
-   if( !env.config() )
-   {
-      WALBERLA_ABORT_NO_DEBUG_INFO( "USAGE: " << argv[0] << " INPUT_FILE" );
-   }
+   Environment env(argc, argv);
+   if (!env.config()) { WALBERLA_ABORT_NO_DEBUG_INFO("USAGE: " << argv[0] << " INPUT_FILE") }
 
    mpi::MPIManager::instance()->useWorldComm();
 
-   const auto & config = *( env.config() );
+   const auto& config = *(env.config());
 
-   Config::BlockHandle configBlock = config.getOneBlock( "ComplexGeometry" );
+   Config::BlockHandle configBlock = config.getOneBlock("ComplexGeometry");
 
-   const std::string     meshFile            = configBlock.getParameter< std::string     >( "meshFile"                                 );
-   const real_t          dx                  = configBlock.getParameter< real_t          >( "coarseDx"                                 );
-   const real_t          omega               = configBlock.getParameter< real_t          >( "coarseOmega"                              );
-   //const uint_t          blockPerProcess     = configBlock.getParameter< uint_t          >( "blocksPerProcess",    uint_t(6)           );
-   const uint_t          timeSteps           = configBlock.getParameter< uint_t          >( "coarseTimeSteps"                          );
-   const Vector3<real_t> bodyForce           = configBlock.getParameter< Vector3<real_t> >( "bodyForce"                                );
-   //const bool            sparseCommunication = configBlock.getParameter< bool            >( "sparseCommunication", true                );
-   const Vector3<real_t> domainBlowUp        = configBlock.getParameter< Vector3<real_t> >( "domainBlowUp",        Vector3<real_t>(6)  );
-   const Vector3<uint_t> blockSize           = configBlock.getParameter< Vector3<uint_t> >( "blockSize",           Vector3<uint_t>(16) );
-         uint_t          numLevels           = configBlock.getParameter< uint_t          >( "numLevels",           uint_t(2)           );
+   const std::string meshFile = configBlock.getParameter< std::string >("meshFile");
+   const real_t dx            = configBlock.getParameter< real_t >("coarseDx");
+   const real_t omega         = configBlock.getParameter< real_t >("coarseOmega");
 
-   numLevels = std::max( numLevels, uint_t(1) );
+   if (configBlock.getParameter< bool >("logLevelDetail"))
+      walberla::logging::Logging::instance()->setLogLevel(walberla::logging::Logging::DETAIL);
 
-   //uint_t numProcesses = uint_c( MPIManager::instance()->numProcesses() );
+   const uint_t timeSteps            = configBlock.getParameter< uint_t >("coarseTimeSteps");
+   const Vector3< real_t > bodyForce = configBlock.getParameter< Vector3< real_t > >("bodyForce");
+   const Vector3< real_t > domainBlowUp =
+      configBlock.getParameter< Vector3< real_t > >("domainBlowUp", Vector3< real_t >(6));
+   const Vector3< uint_t > blockSize =
+      configBlock.getParameter< Vector3< uint_t > >("blockSize", Vector3< uint_t >(16));
+   uint_t numLevels                     = configBlock.getParameter< uint_t >("numLevels", uint_t(2));
+   const bool WriteDistanceOctree = configBlock.getParameter< bool >("WriteDistanceOctree", false);
+   const bool WriteSetupForestAndReturn = configBlock.getParameter< bool >("WriteSetupForestAndReturn", false);
 
-   WALBERLA_LOG_DEVEL_VAR_ON_ROOT( meshFile );
+   numLevels = std::max(numLevels, uint_t(1));
+   const real_t fineDX = dx / real_c(std::pow(2, numLevels));
+
+   // uint_t numProcesses = uint_c( MPIManager::instance()->numProcesses() );
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(meshFile)
 
    auto mesh = make_shared< mesh::TriangleMesh >();
    mesh->request_vertex_colors();
-   WALBERLA_LOG_DEVEL_ON_ROOT( "Loading mesh" );
-   mesh::readAndBroadcast( meshFile, *mesh);
-   vertexToFaceColor( *mesh, mesh::TriangleMesh::Color(255,255,255) );
+   WALBERLA_LOG_DEVEL_ON_ROOT("Loading mesh")
+   mesh::readAndBroadcast(meshFile, *mesh);
+   vertexToFaceColor(*mesh, mesh::TriangleMesh::Color(255, 255, 255));
+
+   WALBERLA_LOG_DEVEL_ON_ROOT("Adding distance info to mesh")
+   auto triDist = make_shared< mesh::TriangleDistance< mesh::TriangleMesh > >(mesh);
+   WALBERLA_LOG_DEVEL_ON_ROOT("Building distance octree")
+   auto distanceOctree = make_shared< mesh::DistanceOctree< mesh::TriangleMesh > >(triDist);
+   WALBERLA_LOG_DEVEL_ON_ROOT("done. Octree has height " << distanceOctree->height())
+
+   // write distance octree to file
+   if (WriteDistanceOctree) {
+      distanceOctree->writeVTKOutput("distanceOctree");
+   }
 
-   WALBERLA_LOG_DEVEL_ON_ROOT( "Adding distance info to mesh" );
-   auto triDist = make_shared< mesh::TriangleDistance<mesh::TriangleMesh> >( mesh );
-   WALBERLA_LOG_DEVEL_ON_ROOT( "Building distance octree" );
-   auto distanceOctree = make_shared< mesh::DistanceOctree<mesh::TriangleMesh> >( triDist );
-   WALBERLA_LOG_DEVEL_ON_ROOT( "done. Octree has height " << distanceOctree->height() );
+   auto aabb = computeAABB(*mesh);
+   aabb.scale(domainBlowUp);
 
-   distanceOctree->writeVTKOutput("distanceOctree");
+   mesh::ComplexGeometryStructuredBlockforestCreator bfc(aabb, Vector3< real_t >(dx));
 
-   auto aabb = computeAABB( *mesh );
-   aabb.scale(  domainBlowUp );
+   bfc.setRootBlockExclusionFunction(mesh::makeExcludeMeshInterior(distanceOctree, dx));
+   bfc.setBlockExclusionFunction(mesh::makeExcludeMeshInteriorRefinement(distanceOctree, fineDX));
 
-   mesh::ComplexGeometryStructuredBlockforestCreator bfc( aabb, Vector3<real_t>( dx ), mesh::makeExcludeMeshInterior( distanceOctree, dx ) );
-   auto meshWorkloadMemory = mesh::makeMeshWorkloadMemory( distanceOctree, dx );
+   auto meshWorkloadMemory = mesh::makeMeshWorkloadMemory(distanceOctree, dx);
    meshWorkloadMemory.setInsideCellWorkload(1);
    meshWorkloadMemory.setOutsideCellWorkload(1);
-   bfc.setWorkloadMemorySUIDAssignmentFunction( meshWorkloadMemory );
-   bfc.setPeriodicity( Vector3<bool>(true) );
-   bfc.setRefinementSelectionFunction( makeRefinementSelection( distanceOctree, numLevels - uint_t(1), dx, dx * real_t(1) ) );
+   bfc.setWorkloadMemorySUIDAssignmentFunction(meshWorkloadMemory);
+   bfc.setPeriodicity(Vector3< bool >(true));
+   bfc.setRefinementSelectionFunction(
+      makeRefinementSelection(distanceOctree, numLevels - uint_t(1), dx, dx * real_t(1)));
 
-   auto structuredBlockforest = bfc.createStructuredBlockForest( blockSize );
+   if (WriteSetupForestAndReturn)
+   {
+      WALBERLA_LOG_INFO_ON_ROOT("Setting up SetupBlockForest")
+      auto setupForest = bfc.createSetupBlockForest(blockSize);
+      WALBERLA_LOG_INFO_ON_ROOT("Writing SetupBlockForest to VTK file")
+      WALBERLA_ROOT_SECTION()
+      {
+         setupForest->writeVTKOutput("SetupBlockForest");
+      }
+      WALBERLA_LOG_INFO_ON_ROOT("Stopping program")
+      return EXIT_SUCCESS;
+   }
+
+   auto structuredBlockforest = bfc.createStructuredBlockForest(blockSize);
 
-   typedef lbm::D3Q19<lbm::collision_model::SRT, false, lbm::force_model::SimpleConstant> LatticeModel_T;
+   typedef lbm::D3Q19< lbm::collision_model::SRT, false, lbm::force_model::SimpleConstant > LatticeModel_T;
 
-   using flag_t = walberla::uint8_t;
-   using FlagField_T = FlagField<flag_t>;
-   using PdfField_T = lbm::PdfField<LatticeModel_T>;
+   using flag_t      = walberla::uint8_t;
+   using FlagField_T = FlagField< flag_t >;
+   using PdfField_T  = lbm::PdfField< LatticeModel_T >;
 
-   LatticeModel_T latticeModel{ lbm::collision_model::SRT( omega ), lbm::force_model::SimpleConstant( bodyForce ) };
+   LatticeModel_T latticeModel{ lbm::collision_model::SRT(omega), lbm::force_model::SimpleConstant(bodyForce) };
 
    static const uint_t NUM_GHOSTLAYERS = 4;
 
-   BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( structuredBlockforest, "pdf field", latticeModel, Vector3<real_t>(0), real_t(1), NUM_GHOSTLAYERS, field::fzyx );
-   BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( structuredBlockforest, "flag field", NUM_GHOSTLAYERS );
+   BlockDataID pdfFieldId = lbm::addPdfFieldToStorage(structuredBlockforest, "pdf field", latticeModel,
+                                                      Vector3< real_t >(0), real_t(1), NUM_GHOSTLAYERS, field::fzyx);
+   BlockDataID flagFieldId =
+      field::addFlagFieldToStorage< FlagField_T >(structuredBlockforest, "flag field", NUM_GHOSTLAYERS);
 
-   const FlagUID fluidFlagUID( "Fluid" );
+   const FlagUID fluidFlagUID("Fluid");
    typedef lbm::DefaultBoundaryHandlingFactory< LatticeModel_T, FlagField_T > BHFactory;
 
-   auto boundariesConfig = configBlock.getOneBlock( "Boundaries" );
+   auto boundariesConfig = configBlock.getOneBlock("Boundaries");
 
-   BlockDataID boundaryHandlingId = BHFactory::addBoundaryHandlingToStorage( structuredBlockforest, "boundary handling", flagFieldId, pdfFieldId, fluidFlagUID,
-                                                                             boundariesConfig.getParameter< Vector3<real_t> >( "velocity0", Vector3<real_t>() ),
-                                                                             boundariesConfig.getParameter< Vector3<real_t> >( "velocity1", Vector3<real_t>() ),
-                                                                             boundariesConfig.getParameter< real_t > ( "pressure0", real_c( 1.0 ) ),
-                                                                             boundariesConfig.getParameter< real_t > ( "pressure1", real_c( 1.001 ) ) );
+   BlockDataID boundaryHandlingId = BHFactory::addBoundaryHandlingToStorage(
+      structuredBlockforest, "boundary handling", flagFieldId, pdfFieldId, fluidFlagUID,
+      boundariesConfig.getParameter< Vector3< real_t > >("velocity0", Vector3< real_t >()),
+      boundariesConfig.getParameter< Vector3< real_t > >("velocity1", Vector3< real_t >()),
+      boundariesConfig.getParameter< real_t >("pressure0", real_c(1.0)),
+      boundariesConfig.getParameter< real_t >("pressure1", real_c(1.001)));
 
-   mesh::ColorToBoundaryMapper<mesh::TriangleMesh> colorToBoundryMapper(( mesh::BoundaryInfo( BHFactory::getNoSlipBoundaryUID() ) ));
+   mesh::ColorToBoundaryMapper< mesh::TriangleMesh > colorToBoundryMapper(
+      (mesh::BoundaryInfo(BHFactory::getNoSlipBoundaryUID())));
 
-  // colorToBoundryMapper.set( mesh::TriangleMesh::Color(255,0,0), mesh::BoundaryInfo( BHFactory::getPressure0BoundaryUID() ) );
-  // colorToBoundryMapper.set( mesh::TriangleMesh::Color(0,0,255), mesh::BoundaryInfo( BHFactory::getPressure1BoundaryUID() ) );
-  // colorToBoundryMapper.set( mesh::TriangleMesh::Color(255,255,255), mesh::BoundaryInfo( BHFactory::getNoSlipBoundaryUID() ) );
+   // colorToBoundryMapper.set( mesh::TriangleMesh::Color(255,0,0), mesh::BoundaryInfo(
+   // BHFactory::getPressure0BoundaryUID() ) ); colorToBoundryMapper.set( mesh::TriangleMesh::Color(0,0,255),
+   // mesh::BoundaryInfo( BHFactory::getPressure1BoundaryUID() ) ); colorToBoundryMapper.set(
+   // mesh::TriangleMesh::Color(255,255,255), mesh::BoundaryInfo( BHFactory::getNoSlipBoundaryUID() ) );
 
-   auto boundaryLocations = colorToBoundryMapper.addBoundaryInfoToMesh( *mesh );
+   auto boundaryLocations = colorToBoundryMapper.addBoundaryInfoToMesh(*mesh);
 
-   mesh::VTKMeshWriter< mesh::TriangleMesh > meshWriter( mesh, "meshBoundaries", 1 );
-   meshWriter.addDataSource( make_shared< mesh::BoundaryUIDFaceDataSource< mesh::TriangleMesh > >( boundaryLocations ) );
-   meshWriter.addDataSource( make_shared< mesh::ColorFaceDataSource< mesh::TriangleMesh > >() );
-   meshWriter.addDataSource( make_shared< mesh::ColorVertexDataSource< mesh::TriangleMesh > >() );
+   mesh::VTKMeshWriter< mesh::TriangleMesh > meshWriter(mesh, "meshBoundaries", 1);
+   meshWriter.addDataSource(make_shared< mesh::BoundaryUIDFaceDataSource< mesh::TriangleMesh > >(boundaryLocations));
+   meshWriter.addDataSource(make_shared< mesh::ColorFaceDataSource< mesh::TriangleMesh > >());
+   meshWriter.addDataSource(make_shared< mesh::ColorVertexDataSource< mesh::TriangleMesh > >());
    meshWriter();
 
-   WALBERLA_LOG_DEVEL_ON_ROOT( "Voxelizing mesh" );
-   mesh::BoundarySetup boundarySetup( structuredBlockforest, makeMeshDistanceFunction( distanceOctree ), NUM_GHOSTLAYERS );
-   //WALBERLA_LOG_DEVEL( "Writing Voxelisation" );
-   //boundarySetup.writeVTKVoxelfile();
-   WALBERLA_LOG_DEVEL_ON_ROOT( "Setting up fluid cells" );
-   boundarySetup.setDomainCells<BHFactory::BoundaryHandling>( boundaryHandlingId, mesh::BoundarySetup::OUTSIDE );
-   WALBERLA_LOG_DEVEL_ON_ROOT( "Setting up boundaries" );
-   boundarySetup.setBoundaries<BHFactory::BoundaryHandling>( boundaryHandlingId, makeBoundaryLocationFunction( distanceOctree, boundaryLocations ), mesh::BoundarySetup::INSIDE );
-   WALBERLA_LOG_DEVEL_ON_ROOT( "done" );
+   WALBERLA_LOG_DEVEL_ON_ROOT("Voxelizing mesh")
+   mesh::BoundarySetup boundarySetup(structuredBlockforest, makeMeshDistanceFunction(distanceOctree), NUM_GHOSTLAYERS);
+   // WALBERLA_LOG_DEVEL( "Writing Voxelisation" );
+   // boundarySetup.writeVTKVoxelfile();
+   WALBERLA_LOG_DEVEL_ON_ROOT("Setting up fluid cells")
+   boundarySetup.setDomainCells< BHFactory::BoundaryHandling >(boundaryHandlingId, mesh::BoundarySetup::OUTSIDE);
+   WALBERLA_LOG_DEVEL_ON_ROOT("Setting up boundaries")
+   boundarySetup.setBoundaries< BHFactory::BoundaryHandling >(
+      boundaryHandlingId, makeBoundaryLocationFunction(distanceOctree, boundaryLocations), mesh::BoundarySetup::INSIDE);
+   WALBERLA_LOG_DEVEL_ON_ROOT("done")
 
-   lbm::BlockForestEvaluation<FlagField_T>( structuredBlockforest, flagFieldId, fluidFlagUID ).logInfoOnRoot();
-   lbm::PerformanceLogger<FlagField_T> perfLogger( structuredBlockforest, flagFieldId, fluidFlagUID, 100 );
+   lbm::BlockForestEvaluation< FlagField_T >(structuredBlockforest, flagFieldId, fluidFlagUID).logInfoOnRoot();
+   lbm::PerformanceLogger< FlagField_T > perfLogger(structuredBlockforest, flagFieldId, fluidFlagUID, 100);
 
-   SweepTimeloop timeloop( structuredBlockforest->getBlockStorage(), timeSteps );
+   SweepTimeloop timeloop(structuredBlockforest->getBlockStorage(), timeSteps);
 
-   auto sweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldId, flagFieldId, fluidFlagUID );
-   auto refinementTimeStep = lbm::refinement::makeTimeStep<LatticeModel_T, BHFactory::BoundaryHandling > ( structuredBlockforest, sweep, pdfFieldId, boundaryHandlingId );
-   timeloop.addFuncBeforeTimeStep( makeSharedFunctor( refinementTimeStep ), "Refinement time step" );   
+   auto sweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >(pdfFieldId, flagFieldId, fluidFlagUID);
+   auto refinementTimeStep = lbm::refinement::makeTimeStep< LatticeModel_T, BHFactory::BoundaryHandling >(
+      structuredBlockforest, sweep, pdfFieldId, boundaryHandlingId);
+   timeloop.addFuncBeforeTimeStep(makeSharedFunctor(refinementTimeStep), "Refinement time step");
 
    // log remaining time
-   timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "remaining time logger" );
+   timeloop.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps()), "remaining time logger");
 
    // LBM stability check
-   timeloop.addFuncAfterTimeStep( makeSharedFunctor( field::makeStabilityChecker< PdfField_T, FlagField_T >( env.config(), structuredBlockforest, pdfFieldId,
-                                                                                                             flagFieldId, fluidFlagUID ) ),
-                                  "LBM stability check" );
-
-   timeloop.addFuncAfterTimeStep( perfLogger, "PerformanceLogger" );
+   timeloop.addFuncAfterTimeStep(makeSharedFunctor(field::makeStabilityChecker< PdfField_T, FlagField_T >(
+                                    env.config(), structuredBlockforest, pdfFieldId, flagFieldId, fluidFlagUID)),
+                                 "LBM stability check");
 
+   timeloop.addFuncAfterTimeStep(perfLogger, "Evaluator: performance logging");
 
    // add VTK output to time loop
-   lbm::VTKOutput< LatticeModel_T, FlagField_T >::addToTimeloop( timeloop, structuredBlockforest, env.config(), pdfFieldId, flagFieldId, fluidFlagUID );
+   lbm::VTKOutput< LatticeModel_T, FlagField_T >::addToTimeloop(timeloop, structuredBlockforest, env.config(),
+                                                                pdfFieldId, flagFieldId, fluidFlagUID);
 
    WcTimingPool timingPool;
-   WALBERLA_LOG_INFO_ON_ROOT( "Starting timeloop" );
-   timeloop.run( timingPool );
-   WALBERLA_LOG_INFO_ON_ROOT( "Timeloop done" );
+   WALBERLA_LOG_INFO_ON_ROOT("Starting timeloop")
+   timeloop.run(timingPool);
+   WALBERLA_LOG_INFO_ON_ROOT("Timeloop done")
    timingPool.unifyRegisteredTimersAcrossProcesses();
-   timingPool.logResultOnRoot( timing::REDUCE_TOTAL, true );
+   timingPool.logResultOnRoot(timing::REDUCE_TOTAL, true);
 
    return EXIT_SUCCESS;
 }
 
-
 } // namespace walberla
 
-int main( int argc, char * argv[] )
-{
-   return walberla::main( argc, argv );
-}
+int main(int argc, char* argv[]) { return walberla::main(argc, argv); }

apps/benchmarks/ComplexGeometry/test.conf

 ComplexGeometry
 {
-   meshFile        cube.obj;
-   coarseDx        0.1;
+   meshFile        bunny.obj;
+   coarseDx        4;
    coarseOmega     1.6;
    coarseTimeSteps 1;
-   numLevels       2;
+   numLevels       4;
    bodyForce       <0.0001, 0, 0>;
    blockSize       <16,16,16>;
    domainBlowUp    <5,5,5>; // simulation domain is blow up factor times mesh size per dimension
+
+   WriteDistanceOctree false;
+   WriteSetupForestAndReturn true;
+   logLevelDetail false;
    
    Boundaries {
    }

apps/benchmarks/CouetteFlow/CMakeLists.txt

 
 waLBerla_link_files_to_builddir( "*.dat" )                 
-                                  
-waLBerla_add_executable( NAME CouetteFlow DEPENDS blockforest boundary core field lbm postprocessing stencil timeloop vtk sqlite )
+
+waLBerla_add_executable( NAME CouetteFlow DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::field walberla::lbm walberla::postprocessing walberla::stencil walberla::timeloop walberla::vtk walberla::sqlite )
 
 ##############
 # Some tests #
 ##############
 
-waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTestNoCheckRelease COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/TestNoCheck.dat --trt --linear-exp PROCESSES 4 CONFIGURATIONS Release RelWithDbgInfo )
-waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTestNoCheckDebug   COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/TestNoCheck.dat --trt --linear-exp PROCESSES 4 LABELS longrun CONFIGURATIONS Debug DebugOptimized )
+waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTestNoCheckRelease COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/TestNoCheck.dat --trt --linear-exp PROCESSES 4 CONFIGURATIONS Release RelWithDbgInfo DEPENDS_ON_TARGETS CouetteFlow )
+waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTestNoCheckDebug   COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/TestNoCheck.dat --trt --linear-exp PROCESSES 4 LABELS longrun CONFIGURATIONS Debug DebugOptimized DEPENDS_ON_TARGETS CouetteFlow )
 
-waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTest0 COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/Test0.dat --trt --linear-exp LABELS longrun CONFIGURATIONS Release RelWithDbgInfo )
-waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTest2 COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/Test2.dat --trt --linear-exp LABELS longrun verylongrun PROCESSES 4 CONFIGURATIONS Release RelWithDbgInfo )
+waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTest0 COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/Test0.dat --trt --linear-exp LABELS longrun CONFIGURATIONS Release RelWithDbgInfo DEPENDS_ON_TARGETS CouetteFlow )
+waLBerla_execute_test( NO_MODULE_LABEL NAME CouetteFlowTest2 COMMAND $<TARGET_FILE:CouetteFlow> ${CMAKE_CURRENT_SOURCE_DIR}/Test2.dat --trt --linear-exp LABELS longrun verylongrun PROCESSES 4 CONFIGURATIONS Release RelWithDbgInfo DEPENDS_ON_TARGETS CouetteFlow )

apps/benchmarks/CouetteFlow/CouetteFlow.cpp

@@ -162,19 +162,19 @@ template< typename LatticeModel_T, class Enable = void >
 struct StencilString;
 
 template< typename LatticeModel_T >
-struct StencilString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::Stencil, stencil::D3Q15 >::value >::type >
+struct StencilString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::Stencil, stencil::D3Q15 > > >
 {
    static const char * str() { return "D3Q15"; }
 };
 
 template< typename LatticeModel_T >
-struct StencilString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::Stencil, stencil::D3Q19 >::value >::type >
+struct StencilString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::Stencil, stencil::D3Q19 > > >
 {
    static const char * str() { return "D3Q19"; }
 };
 
 template< typename LatticeModel_T >
-struct StencilString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::Stencil, stencil::D3Q27 >::value >::type >
+struct StencilString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::Stencil, stencil::D3Q27 > > >
 {
    static const char * str() { return "D3Q27"; }
 };
@@ -184,22 +184,22 @@ template< typename LatticeModel_T, class Enable = void >
 struct CollisionModelString;
 
 template< typename LatticeModel_T >
-struct CollisionModelString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::CollisionModel::tag,
-                                                                                          lbm::collision_model::SRT_tag >::value >::type >
+struct CollisionModelString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::CollisionModel::tag,
+                                                                                          lbm::collision_model::SRT_tag > > >
 {
    static const char * str() { return "SRT"; }
 };
 
 template< typename LatticeModel_T >
-struct CollisionModelString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::CollisionModel::tag,
-                                                                                          lbm::collision_model::TRT_tag >::value >::type >
+struct CollisionModelString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::CollisionModel::tag,
+                                                                                          lbm::collision_model::TRT_tag > > >
 {
    static const char * str() { return "TRT"; }
 };
 
 template< typename LatticeModel_T >
-struct CollisionModelString< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::CollisionModel::tag,
-                                                                                          lbm::collision_model::MRT_tag >::value >::type >
+struct CollisionModelString< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::CollisionModel::tag,
+                                                                                          lbm::collision_model::MRT_tag > > >
 {
    static const char * str() { return "MRT"; }
 };
@@ -271,10 +271,10 @@ private:
 
 static void workloadAndMemoryAssignment( SetupBlockForest& forest, const memory_t memoryPerBlock )
 {
-   for( auto block = forest.begin(); block != forest.end(); ++block )
+   for(auto & block : forest)
    {
-      block->setWorkload( numeric_cast< workload_t >( uint_t(1) << block->getLevel() ) );
-      block->setMemory( memoryPerBlock );
+      block.setWorkload( numeric_cast< workload_t >( uint_t(1) << block.getLevel() ) );
+      block.setMemory( memoryPerBlock );
    }
 }
 
@@ -293,7 +293,7 @@ static shared_ptr< SetupBlockForest > createSetupBlockForest( const blockforest:
                                                              ( setup.zCells + uint_t(2) * FieldGhostLayers ) ) * memoryPerCell;
 
    forest->addRefinementSelectionFunction( refinementSelectionFunctions );
-   forest->addWorkloadMemorySUIDAssignmentFunction( std::bind( workloadAndMemoryAssignment, std::placeholders::_1, memoryPerBlock ) );
+   forest->addWorkloadMemorySUIDAssignmentFunction([memoryPerBlock](auto & sbf) { workloadAndMemoryAssignment(sbf, memoryPerBlock); });
 
    forest->init( AABB( real_c(0), real_c(0), real_c(0), real_c( setup.xBlocks * setup.xCells ),
                                                         real_c( setup.yBlocks * setup.yCells ),
@@ -634,10 +634,10 @@ struct AddRefinementTimeStep
 };
 
 template< typename LatticeModel_T  >
-struct AddRefinementTimeStep< LatticeModel_T, typename std::enable_if< std::is_same< typename LatticeModel_T::CollisionModel::tag, lbm::collision_model::MRT_tag >::value ||
-                                                                       std::is_same< typename LatticeModel_T::Stencil, stencil::D3Q15 >::value ||
-                                                                       std::is_same< typename LatticeModel_T::Stencil, stencil::D3Q27 >::value
-                                                                       >::type >
+struct AddRefinementTimeStep< LatticeModel_T, std::enable_if_t< std::is_same_v< typename LatticeModel_T::CollisionModel::tag, lbm::collision_model::MRT_tag > ||
+                                                                       std::is_same_v< typename LatticeModel_T::Stencil, stencil::D3Q15 > ||
+                                                                       std::is_same_v< typename LatticeModel_T::Stencil, stencil::D3Q27 >
+                                                                       > >
 {
    static void add( SweepTimeloop & timeloop, shared_ptr< blockforest::StructuredBlockForest > & blocks,
                     const BlockDataID & pdfFieldId, const BlockDataID & flagFieldId, const BlockDataID & boundaryHandlingId,
@@ -718,9 +718,9 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
 
    if( vtkBeforeTimeStep )
    {
-      for( auto output = vtkOutputFunctions.begin(); output != vtkOutputFunctions.end(); ++output )
-         timeloop.addFuncBeforeTimeStep( output->second.outputFunction, std::string("VTK: ") + output->first,
-                                         output->second.requiredGlobalStates, output->second.incompatibleGlobalStates );
+      for(auto & output : vtkOutputFunctions)
+         timeloop.addFuncBeforeTimeStep( output.second.outputFunction, std::string("VTK: ") + output.first,
+                                        output.second.requiredGlobalStates, output.second.incompatibleGlobalStates );
    }
 
    // add 'refinement' LB time step to time loop
@@ -733,11 +733,11 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
 
    // evaluation
 
-   const auto exactSolutionFunction = std::bind( exactVelocity, std::placeholders::_1, blocks->getDomain(), setup.maxVelocity_L );
+   const auto exactSolutionFunction = [domain = blocks->getDomain(), maxVel = setup.maxVelocity_L](auto & p) { return exactVelocity(p, domain, maxVel); };
 
    auto volumetricFlowRate = field::makeVolumetricFlowRateEvaluation< VelocityAdaptor_T, FlagField_T >( configBlock, blocks, velocityAdaptorId,
                                                                                                         flagFieldId, Fluid_Flag,
-                                                                                                        std::bind( exactFlowRate, setup.flowRate_L ),
+                                                                                                        [flowRate = setup.flowRate_L] { return exactFlowRate(flowRate); },
                                                                                                         exactSolutionFunction );
    volumetricFlowRate->setNormalizationFactor( real_t(1) / setup.maxVelocity_L );
    volumetricFlowRate->setDomainNormalization( Vector3<real_t>( real_t(1) ) );
@@ -766,14 +766,14 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
 
    if( !vtkBeforeTimeStep )
    {
-      for( auto output = vtkOutputFunctions.begin(); output != vtkOutputFunctions.end(); ++output )
-         timeloop.addFuncAfterTimeStep( output->second.outputFunction, std::string("VTK: ") + output->first,
-                                        output->second.requiredGlobalStates, output->second.incompatibleGlobalStates );
+      for(auto & output : vtkOutputFunctions)
+         timeloop.addFuncAfterTimeStep( output.second.outputFunction, std::string("VTK: ") + output.first,
+                                       output.second.requiredGlobalStates, output.second.incompatibleGlobalStates );
    }
 
    // remaining time logger
 
-   const double remainingTimeLoggerFrequency = configBlock.getParameter< double >( "remainingTimeLoggerFrequency", 3.0 );
+   const real_t remainingTimeLoggerFrequency = configBlock.getParameter< real_t >( "remainingTimeLoggerFrequency", real_c(3.0) );
    timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger( timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency ), "Remaining time logger" );
 
    // logging right before the simulation starts

apps/benchmarks/DEM/CMakeLists.txt

-waLBerla_add_executable( NAME DEM FILES DEM.cpp DEPENDS blockforest core pe )
+waLBerla_add_executable( NAME DEM FILES DEM.cpp DEPENDS walberla::blockforest walberla::core walberla::pe )

apps/benchmarks/FieldCommunication/CMakeLists.txt

@@ -4,4 +4,4 @@ waLBerla_link_files_to_builddir( "*.py" )
 
 
 waLBerla_add_executable ( NAME FieldCommunication
-                          DEPENDS blockforest core domain_decomposition field postprocessing sqlite python_coupling )
+      DEPENDS walberla::blockforest walberla::core walberla::domain_decomposition walberla::field walberla::postprocessing walberla::sqlite walberla::python_coupling )

apps/benchmarks/FieldCommunication/FieldCommunication.cpp

@@ -34,7 +34,7 @@ template<typename Stencil_T>
 class SingleMessageBufferedScheme
 {
 public:
-    typedef Stencil_T Stencil;
+    using Stencil = Stencil_T;
 
     SingleMessageBufferedScheme( const weak_ptr< StructuredBlockForest > & bf, const int tag = 17953 )
             : blockForest_( bf ), tag_( tag ) {}
@@ -354,8 +354,8 @@ int main( int argc, char **argv )
             auto databaseBlock = config->getBlock( "Database" );
             if ( databaseBlock )
             {
-                for ( auto it = databaseBlock.begin(); it != databaseBlock.end(); ++it )
-                    stringProperties[it->first] = it->second;
+                for (const auto & it : databaseBlock)
+                    stringProperties[it.first] = it.second;
             }
 
             realProperties["total_min"] = real_c( timingPool["totalTime"].min()) / real_c( iterations );

apps/benchmarks/FluidParticleCoupling/CMakeLists.txt

@@ -11,63 +11,61 @@ if( WALBERLA_BUILD_WITH_CODEGEN )
             )
 
    waLBerla_add_executable(NAME SphereWallCollision FILES SphereWallCollision.cpp
-         DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling
-         mesa_pd postprocessing timeloop vtk FluidParticleCouplingGeneratedLBM)
+         DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling
+         walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
 
    waLBerla_add_executable(NAME SettlingSphereInBox FILES SettlingSphereInBox.cpp
-         DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling
-         mesa_pd postprocessing timeloop vtk FluidParticleCouplingGeneratedLBM)
+         DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling
+         walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
 
    waLBerla_add_executable(NAME SphereMovingWithPrescribedVelocity FILES SphereMovingWithPrescribedVelocity.cpp
-         DEPENDS blockforest boundary core domain_decomposition field lbm mesa_pd lbm_mesapd_coupling
-         postprocessing timeloop vtk FluidParticleCouplingGeneratedLBM)
+         DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::mesa_pd walberla::lbm_mesapd_coupling
+         walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
 
    waLBerla_add_executable(NAME LubricationForceEvaluation FILES LubricationForceEvaluation.cpp
-         DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd
-         postprocessing timeloop vtk FluidParticleCouplingGeneratedLBM)
+         DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd
+         walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
 
    waLBerla_add_executable(NAME DragForceSphere FILES DragForceSphere.cpp
-         DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd
-         postprocessing timeloop vtk FluidParticleCouplingGeneratedLBMWithForce)
+         DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd
+         walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBMWithForce )
 
    waLBerla_add_executable(NAME ForcesOnSphereNearPlane FILES ForcesOnSphereNearPlane.cpp
-         DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd
-         postprocessing timeloop vtk FluidParticleCouplingGeneratedLBM)
+         DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd
+         walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
 
    waLBerla_add_executable(NAME ObliqueWetCollision FILES ObliqueWetCollision.cpp
-         DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd
-         postprocessing timeloop vtk FluidParticleCouplingGeneratedLBM)
+         DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd
+         walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
 
    waLBerla_add_executable(NAME MotionSettlingSphere FILES MotionSettlingSphere.cpp
-         DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd
-         postprocessing timeloop vtk FluidParticleCouplingGeneratedLBM)
+         DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd
+         walberla::postprocessing walberla::timeloop walberla::vtk FluidParticleCouplingGeneratedLBM )
 
 else()
 
 waLBerla_add_executable ( NAME SphereWallCollision FILES SphereWallCollision.cpp
-                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+      DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
 
 waLBerla_add_executable ( NAME SettlingSphereInBox FILES SettlingSphereInBox.cpp
-                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+      DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
 
 waLBerla_add_executable ( NAME SphereMovingWithPrescribedVelocity FILES SphereMovingWithPrescribedVelocity.cpp
-                          DEPENDS blockforest boundary core domain_decomposition field lbm mesa_pd lbm_mesapd_coupling postprocessing timeloop vtk )
+      DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::mesa_pd walberla::lbm_mesapd_coupling walberla::postprocessing walberla::timeloop walberla::vtk )
 
 waLBerla_add_executable ( NAME LubricationForceEvaluation FILES LubricationForceEvaluation.cpp
-                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+      DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
 
 waLBerla_add_executable ( NAME DragForceSphere FILES DragForceSphere.cpp
-                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+      DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
 
 waLBerla_add_executable ( NAME ForcesOnSphereNearPlane FILES ForcesOnSphereNearPlane.cpp
-                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+      DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
 
 waLBerla_add_executable ( NAME ObliqueWetCollision FILES ObliqueWetCollision.cpp
-                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+      DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::domain_decomposition walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
 
 endif()
 
-
-
 waLBerla_add_executable ( NAME ObliqueDryCollision FILES ObliqueDryCollision.cpp
-      DEPENDS blockforest core mesa_pd postprocessing )
\ No newline at end of file
+      DEPENDS walberla::blockforest walberla::core walberla::mesa_pd walberla::postprocessing )

apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp

@@ -216,13 +216,8 @@ public:
 
       WALBERLA_MPI_SECTION()
       {
-         mpi::allReduceInplace( force[0], mpi::SUM );
-         mpi::allReduceInplace( force[1], mpi::SUM );
-         mpi::allReduceInplace( force[2], mpi::SUM );
-
-         mpi::allReduceInplace( torque[0], mpi::SUM );
-         mpi::allReduceInplace( torque[1], mpi::SUM );
-         mpi::allReduceInplace( torque[2], mpi::SUM );
+         mpi::allReduceInplace( force, mpi::SUM );
+         mpi::allReduceInplace( torque, mpi::SUM );
       }
 
       if( fileIO_ )

apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py

 import sympy as sp
 from sympy.core.cache import clear_cache
 import pystencils as ps
-from lbmpy.creationfunctions import create_lb_method_from_existing, create_lb_method
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, Method, Stencil, create_lb_method
 from lbmpy_walberla import generate_lattice_model
 from pystencils_walberla import CodeGeneration
 from pystencils_walberla import get_vectorize_instruction_set
 
 from lbmpy.creationfunctions import create_lb_collision_rule
 from lbmpy.moments import is_even, get_order, MOMENT_SYMBOLS
-from lbmpy.stencils import get_stencil
+from lbmpy.stencils import LBStencil
 
 from collections import OrderedDict
 
 with CodeGeneration() as ctx:
 
     generatedMethod = 'TRTlike'
-    #generatedMethod = 'D3Q27TRTlike'
-    #generatedMethod = 'cumulant'
+    # generatedMethod = 'D3Q27TRTlike'
+    # generatedMethod = 'cumulant'
 
     clear_cache()
 
@@ -28,7 +28,7 @@ with CodeGeneration() as ctx:
         omegaVisc = sp.Symbol('omega_visc')
         omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
         omegaMagic = sp.Symbol('omega_magic')
-        stencil = get_stencil('D3Q19', 'walberla')
+        stencil = LBStencil(Stencil.D3Q19)
 
         x, y, z = MOMENT_SYMBOLS
         one = sp.Rational(1, 1)
@@ -44,33 +44,43 @@ with CodeGeneration() as ctx:
         ]
 
         # relaxation rate for first group of moments (1,x,y,z) is set to zero internally
-        relaxation_rates=[omegaBulk.center_vector, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaVisc, omegaMagic]
+        relaxation_rates = [omegaBulk.center, omegaBulk.center,
+                            omegaMagic, omegaVisc, omegaVisc, omegaMagic]
 
-        methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False,
-                                           nested_moments=moments, relaxation_rates=relaxation_rates)
+        lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, continuous_equilibrium=False,
+                               zero_centered=False, delta_equilibrium=False,
+                               nested_moments=moments, relaxation_rates=relaxation_rates)
 
-        #print(methodWithForce.relaxation_rates)
-        #print(methodWithForce.moment_matrix)
-        collision_rule = create_lb_collision_rule(lb_method=methodWithForce, optimization={'cse_global': True})
+        lbm_opt = LBMOptimisation(cse_global=True)
+
+        methodWithForce = create_lb_method(lbm_config=lbm_config)
+
+        # print(methodWithForce.relaxation_rates)
+        # print(methodWithForce.moment_matrix)
+        collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
         generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx', cpu_vectorize_info=cpu_vectorize_info)
 
- 
     if generatedMethod == 'D3Q27TRTlike':
 
         omegaVisc = sp.Symbol('omega_visc')
         omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
         omegaMagic = sp.Symbol('omega_magic')
-        stencil = get_stencil('D3Q27', 'walberla')
+        stencil = LBStencil(Stencil.D3Q27)
 
         relaxation_rates = [omegaVisc, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaMagic, omegaVisc]
 
-        methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False, weighted=True,
-                                           compressible=False, relaxation_rates=relaxation_rates)
+        lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, maxwellian_moments=False, weighted=True,
+                               compressible=False, relaxation_rates=relaxation_rates)
+
+        lbm_opt = LBMOptimisation(cse_global=True)
+
+        methodWithForce = create_lb_method(lbm_config=lbm_config)
 
-        collision_rule = create_lb_collision_rule(lb_method=methodWithForce, optimization={'cse_global': True})
+        collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
         generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx',
                                cpu_vectorize_info=cpu_vectorize_info)
 
+# using create_with_continuous_maxwellian_eq_moments won't work probably
     if generatedMethod == 'cumulant':
 
         x, y, z = MOMENT_SYMBOLS
@@ -122,7 +132,7 @@ with CodeGeneration() as ctx:
             else:
                 return 1
 
-        stencil = get_stencil('D3Q27', 'walberla')
+        stencil = LBStencil(Stencil.D3Q27)
 
         omega = sp.Symbol('omega')
         rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
@@ -139,6 +149,3 @@ with CodeGeneration() as ctx:
 
         generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx',
                                cpu_vectorize_info=cpu_vectorize_info)
-
-
-

apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py

 import sympy as sp
 from sympy.core.cache import clear_cache
 import pystencils as ps
-from lbmpy.creationfunctions import create_lb_method_from_existing, create_lb_ast, create_lb_method
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, ForceModel, Method, Stencil, create_lb_method
 from lbmpy_walberla import generate_lattice_model
 from pystencils_walberla import CodeGeneration
 from pystencils_walberla import get_vectorize_instruction_set
 
 from lbmpy.creationfunctions import create_lb_collision_rule
 from lbmpy.moments import MOMENT_SYMBOLS, is_even, get_order
-from lbmpy.stencils import get_stencil
-from lbmpy.forcemodels import Luo
+from lbmpy.stencils import LBStencil
 
 from collections import OrderedDict
 
 with CodeGeneration() as ctx:
 
     forcing = (sp.symbols('fx'), 0, 0)
-    forcemodel = Luo(forcing)
 
     generatedMethod = 'TRTlike'
     # generatedMethod = 'D3Q27TRTlike'
@@ -35,7 +33,7 @@ with CodeGeneration() as ctx:
         omegaVisc = sp.Symbol('omega_visc')
         omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
         omegaMagic = sp.Symbol('omega_magic')
-        stencil = get_stencil('D3Q19', 'walberla')
+        stencil = LBStencil(Stencil.D3Q19)
 
         x, y, z = MOMENT_SYMBOLS
         one = sp.Rational(1, 1)
@@ -51,29 +49,42 @@ with CodeGeneration() as ctx:
         ]
 
         # relaxation rate for first group of moments (1,x,y,z) is set to zero internally
-        relaxation_rates=[omegaBulk.center_vector, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaVisc, omegaMagic]
+        relaxation_rates = [omegaBulk.center, omegaBulk.center,
+                            omegaMagic, omegaVisc, omegaVisc, omegaMagic]
 
-        methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False,
-                                           force_model=forcemodel, nested_moments=moments, relaxation_rates=relaxation_rates)
+        lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, continuous_equilibrium=False,
+                               zero_centered=False, delta_equilibrium=False,
+                               force=forcing, force_model=ForceModel.LUO,
+                               nested_moments=moments, relaxation_rates=relaxation_rates)
+
+        lbm_opt = LBMOptimisation(cse_global=True)
+
+        methodWithForce = create_lb_method(lbm_config=lbm_config)
+
+        # print(methodWithForce.relaxation_rates)
+        # print(methodWithForce.moment_matrix)
+        collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+        generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
+                               cpu_vectorize_info=cpu_vectorize_info)
 
-        #print(methodWithForce.relaxation_rates)
-        #print(methodWithForce.moment_matrix)
-        collision_rule = create_lb_collision_rule(lb_method=methodWithForce, optimization={'cse_global': True})
-        generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx', cpu_vectorize_info=cpu_vectorize_info)
- 
     if generatedMethod == 'D3Q27TRTlike':
 
         omegaVisc = sp.Symbol('omega_visc')
         omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
         omegaMagic = sp.Symbol('omega_magic')
-        stencil = get_stencil('D3Q27', 'walberla')
+        stencil = LBStencil(Stencil.D3Q27)
 
         relaxation_rates = [omegaVisc, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaMagic, omegaVisc]
 
-        methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False, weighted=True,
-                                           compressible=False, force_model=forcemodel, relaxation_rates=relaxation_rates)
+        lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, maxwellian_moments=False, weighted=True,
+                               force=forcing, force_model=ForceModel.LUO,
+                               compressible=False, relaxation_rates=relaxation_rates)
+
+        lbm_opt = LBMOptimisation(cse_global=True)
 
-        collision_rule = create_lb_collision_rule(lb_method=methodWithForce, optimization={'cse_global': True})
+        methodWithForce = create_lb_method(lbm_config=lbm_config)
+
+        collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
         generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
                                cpu_vectorize_info=cpu_vectorize_info)
 
@@ -128,7 +139,7 @@ with CodeGeneration() as ctx:
             else:
                 return 1
 
-        stencil = get_stencil('D3Q27', 'walberla')
+        stencil = LBStencil(Stencil.D3Q27)
 
         omega = sp.Symbol('omega')
         rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
@@ -199,7 +210,7 @@ with CodeGeneration() as ctx:
             else:
                 return omegaMagic
 
-        stencil = get_stencil('D3Q27', 'walberla')
+        stencil = LBStencil(Stencil.D3Q27)
 
         omegaVisc = sp.Symbol('omega_visc')
         omegaMagic = sp.Symbol('omega_magic')
@@ -219,4 +230,3 @@ with CodeGeneration() as ctx:
 
         generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
                                cpu_vectorize_info=cpu_vectorize_info)
-

apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp

@@ -705,18 +705,10 @@ int main( int argc, char **argv )
 
          WALBERLA_MPI_SECTION()
          {
-            mpi::allReduceInplace( hydForce[0], mpi::SUM );
-            mpi::allReduceInplace( hydForce[1], mpi::SUM );
-            mpi::allReduceInplace( hydForce[2], mpi::SUM );
-            mpi::reduceInplace( lubForce[0], mpi::SUM );
-            mpi::reduceInplace( lubForce[1], mpi::SUM );
-            mpi::reduceInplace( lubForce[2], mpi::SUM );
-            mpi::allReduceInplace( hydTorque[0], mpi::SUM );
-            mpi::allReduceInplace( hydTorque[1], mpi::SUM );
-            mpi::allReduceInplace( hydTorque[2], mpi::SUM );
-            mpi::reduceInplace( lubTorque[0], mpi::SUM );
-            mpi::reduceInplace( lubTorque[1], mpi::SUM );
-            mpi::reduceInplace( lubTorque[2], mpi::SUM );
+            mpi::allReduceInplace( hydForce, mpi::SUM );
+            mpi::reduceInplace( lubForce, mpi::SUM );
+            mpi::allReduceInplace( hydTorque, mpi::SUM );
+            mpi::reduceInplace( lubTorque, mpi::SUM );
          }
 
          curForceNorm = hydForce.length();

apps/benchmarks/FluidParticleCoupling/MotionSettlingSphere.cpp

@@ -313,9 +313,7 @@ private:
 
       WALBERLA_MPI_SECTION()
       {
-         mpi::allReduceInplace( force[0], mpi::SUM );
-         mpi::allReduceInplace( force[1], mpi::SUM );
-         mpi::allReduceInplace( force[2], mpi::SUM );
+         mpi::allReduceInplace( force, mpi::SUM );
       }
       return force;
    }
@@ -393,41 +391,26 @@ public:
    // evaluate and write the sphere properties
    void operator()(const uint_t timestep)
    {
-      Vector3<real_t> transVel( real_t(0) );
-      Vector3<real_t> angularVel( real_t(0) );
-      Vector3<real_t> pos( real_t(0) );
+      Vector3<real_t> particleTransVel( real_t(0) );
+      Vector3<real_t> particleAngularVel( real_t(0) );
+      Vector3<real_t> particlePos( real_t(0) );
 
-      Vector3<real_t> force( real_t(0) );
-      Vector3<real_t> torque( real_t(0) );
+      Vector3<real_t> particleForce( real_t(0) );
+      Vector3<real_t> particleTorque( real_t(0) );
 
       size_t idx = ac_->uidToIdx(sphereUid_);
       if( idx != ac_->getInvalidIdx())
       {
          if(!mesa_pd::data::particle_flags::isSet( ac_->getFlags(idx), mesa_pd::data::particle_flags::GHOST))
          {
-            pos = ac_->getPosition(idx);
-            transVel = ac_->getLinearVelocity(idx);
-            angularVel = ac_->getAngularVelocity(idx);
-            force = ac_->getHydrodynamicForce(idx);
-            torque = ac_->getHydrodynamicTorque(idx);
+            particlePos = ac_->getPosition(idx);
+            particleTransVel = ac_->getLinearVelocity(idx);
+            particleAngularVel = ac_->getAngularVelocity(idx);
+            particleForce = ac_->getHydrodynamicForce(idx);
+            particleTorque = ac_->getHydrodynamicTorque(idx);
          }
       }
 
-      std::vector<real_t> particlePos(3);
-      particlePos[0]=pos[0]; particlePos[1]=pos[1]; particlePos[2]=pos[2];
-
-      std::vector<real_t> particleTransVel(3);
-      particleTransVel[0]=transVel[0]; particleTransVel[1]=transVel[1]; particleTransVel[2]=transVel[2];
-
-      std::vector<real_t> particleAngularVel(3);
-      particleAngularVel[0]=angularVel[0]; particleAngularVel[1]=angularVel[1]; particleAngularVel[2]=angularVel[2];
-
-      std::vector<real_t> particleForce(3);
-      particleForce[0]=force[0]; particleForce[1]=force[1]; particleForce[2]=force[2];
-
-      std::vector<real_t> particleTorque(3);
-      particleTorque[0]=torque[0]; particleTorque[1]=torque[1]; particleTorque[2]=torque[2];
-
       // reduce to root
       WALBERLA_MPI_SECTION()
       {
@@ -547,10 +530,7 @@ public:
 
       WALBERLA_MPI_SECTION()
       {
-         mpi::allReduceInplace( u_p[0], mpi::SUM );
-         mpi::allReduceInplace( u_p[1], mpi::SUM );
-         mpi::allReduceInplace( u_p[2], mpi::SUM );
-
+         mpi::allReduceInplace( u_p, mpi::SUM );
       }
 
       Vector3<real_t> u_p_r = u_p - u_infty_;

apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp

@@ -268,15 +268,9 @@ public:
 
       WALBERLA_MPI_SECTION()
       {
-         mpi::allReduceInplace( pos[0], mpi::SUM );
-         mpi::allReduceInplace( pos[1], mpi::SUM );
-         mpi::allReduceInplace( pos[2], mpi::SUM );
-         mpi::allReduceInplace( transVel[0], mpi::SUM );
-         mpi::allReduceInplace( transVel[1], mpi::SUM );
-         mpi::allReduceInplace( transVel[2], mpi::SUM );
-         mpi::allReduceInplace( angVel[0], mpi::SUM );
-         mpi::allReduceInplace( angVel[1], mpi::SUM );
-         mpi::allReduceInplace( angVel[2], mpi::SUM );
+         mpi::allReduceInplace( pos, mpi::SUM );
+         mpi::allReduceInplace( transVel, mpi::SUM );
+         mpi::allReduceInplace( angVel, mpi::SUM );
       }
 
       position_ = pos;
@@ -451,9 +445,7 @@ Vector3<real_t> getForce(walberla::id_t uid, ParticleAccessor_T & ac)
    }
    WALBERLA_MPI_SECTION()
    {
-      mpi::allReduceInplace( force[0], mpi::SUM );
-      mpi::allReduceInplace( force[1], mpi::SUM );
-      mpi::allReduceInplace( force[2], mpi::SUM );
+      mpi::allReduceInplace( force, mpi::SUM );
    }
    return force;
 }
@@ -469,9 +461,7 @@ Vector3<real_t> getTorque(walberla::id_t uid, ParticleAccessor_T & ac)
    }
    WALBERLA_MPI_SECTION()
    {
-      mpi::allReduceInplace( torque[0], mpi::SUM );
-      mpi::allReduceInplace( torque[1], mpi::SUM );
-      mpi::allReduceInplace( torque[2], mpi::SUM );
+      mpi::allReduceInplace( torque, mpi::SUM );
    }
    return torque;
 }

apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp

@@ -223,17 +223,9 @@ public:
 
       WALBERLA_MPI_SECTION()
       {
-         mpi::allReduceInplace( pos[0], mpi::SUM );
-         mpi::allReduceInplace( pos[1], mpi::SUM );
-         mpi::allReduceInplace( pos[2], mpi::SUM );
-
-         mpi::allReduceInplace( transVel[0], mpi::SUM );
-         mpi::allReduceInplace( transVel[1], mpi::SUM );
-         mpi::allReduceInplace( transVel[2], mpi::SUM );
-
-         mpi::allReduceInplace( hydForce[0], mpi::SUM );
-         mpi::allReduceInplace( hydForce[1], mpi::SUM );
-         mpi::allReduceInplace( hydForce[2], mpi::SUM );
+         mpi::allReduceInplace( pos, mpi::SUM );
+         mpi::allReduceInplace( transVel, mpi::SUM );
+         mpi::allReduceInplace( hydForce, mpi::SUM );
       }
 
       position_ = pos[2];

apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp

@@ -361,9 +361,7 @@ Vector3<real_t> getVelocityAtPosition(const shared_ptr<StructuredBlockStorage> &
       }
 
    }
-   mpi::reduceInplace(vel[0], mpi::SUM);
-   mpi::reduceInplace(vel[1], mpi::SUM);
-   mpi::reduceInplace(vel[2], mpi::SUM);
+   mpi::reduceInplace(vel, mpi::SUM);
 
    return vel;
 }
@@ -878,7 +876,7 @@ int main( int argc, char **argv )
       real_t defaultOmegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, real_t(1));
       shared_ptr<OmegaBulkAdapter_T> omegaBulkAdapter = make_shared<OmegaBulkAdapter_T>(blocks, omegaBulkFieldID, accessor, defaultOmegaBulk, omegaBulk, adaptionLayerSize, sphereSelector);
       timeloopAfterParticles.add() << Sweep( makeSharedSweep(omegaBulkAdapter), "Omega Bulk Adapter");
-      // initally adapt
+      // initially adapt
       for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
          (*omegaBulkAdapter)(blockIt.get());
       }

apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp

@@ -264,9 +264,7 @@ public:
 
       WALBERLA_MPI_SECTION()
       {
-         mpi::allReduceInplace( pos[0], mpi::SUM );
-         mpi::allReduceInplace( pos[1], mpi::SUM );
-         mpi::allReduceInplace( pos[2], mpi::SUM );
+         mpi::allReduceInplace( pos, mpi::SUM );
          mpi::allReduceInplace( transVel[2], mpi::SUM );
       }
 

apps/benchmarks/FluidParticleCouplingWithLoadBalancing/CMakeLists.txt

-waLBerla_add_executable( NAME FluidParticleWorkloadEvaluation FILES FluidParticleWorkloadEvaluation.cpp DEPENDS blockforest boundary core field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+waLBerla_add_executable( NAME FluidParticleWorkloadEvaluation FILES FluidParticleWorkloadEvaluation.cpp DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )
 
-waLBerla_add_executable( NAME FluidParticleWorkloadDistribution FILES FluidParticleWorkloadDistribution.cpp DEPENDS blockforest boundary core field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+waLBerla_add_executable( NAME FluidParticleWorkloadDistribution FILES FluidParticleWorkloadDistribution.cpp DEPENDS walberla::blockforest walberla::boundary walberla::core walberla::field walberla::lbm walberla::lbm_mesapd_coupling walberla::mesa_pd walberla::postprocessing walberla::timeloop walberla::vtk )

apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadDistribution.cpp

@@ -439,14 +439,14 @@ void evaluateParticleSimulation(const shared_ptr<mesa_pd::data::ParticleStorage>
 {
    auto numShapes = ss->shapes.size();
    uint_t numLocalParticles = 0;
-   uint_t numGhostParticles = 0;
+   //uint_t numGhostParticles = 0;
    uint_t numGlobalParticlesOfRank = 0;
    for(auto p = ps->begin(); p != ps->end(); ++p)
    {
       using namespace walberla::mesa_pd::data::particle_flags;
       if (isSet(p->getFlags(), GHOST))
       {
-         ++numGhostParticles;
+         //++numGhostParticles;
       } else if (isSet(p->getFlags(), GLOBAL))
       {
          ++numGlobalParticlesOfRank;
@@ -843,7 +843,7 @@ int main( int argc, char **argv )
    auto sphereShape = ss->create<mesa_pd::data::Sphere>( diameter * real_t(0.5) );
    ss->shapes[sphereShape]->updateMassAndInertia(densityRatio);
 
-   std::mt19937 randomGenerator (static_cast<unsigned int>(2610)); // fixed seed: quasi-random and reproducable
+   std::mt19937 randomGenerator (static_cast<unsigned int>(2610)); // fixed seed: quasi-random and reproducible
 
    for( uint_t nSed = 0; nSed < numberOfSediments; ++nSed )
    {
@@ -962,7 +962,7 @@ int main( int argc, char **argv )
 
       if(currentPhase == 1)
       {
-         // damp velocites to avoid too large ones
+         // damp velocities to avoid too large ones
          ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
                               [](const size_t idx, ParticleAccessor_T& ac){
                                  ac.setLinearVelocity(idx, ac.getLinearVelocity(idx) * real_t(0.5));
@@ -1051,14 +1051,14 @@ int main( int argc, char **argv )
    LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega, lbm::collision_model::TRT::threeSixteenth ) );
 
    // add PDF field
-   BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
+   BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (fzyx)", latticeModel,
                                                                          Vector3< real_t >( real_t(0) ), real_t(1),
-                                                                         FieldGhostLayers, field::zyxf );
+                                                                         FieldGhostLayers, field::fzyx );
    // add flag field
    BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field", FieldGhostLayers );
 
    // add particle field
-   BlockDataID particleFieldID = field::addToStorage<lbm_mesapd_coupling::ParticleField_T>( blocks, "particle field", accessor->getInvalidUid(), field::zyxf, FieldGhostLayers );
+   BlockDataID particleFieldID = field::addToStorage<lbm_mesapd_coupling::ParticleField_T>( blocks, "particle field", accessor->getInvalidUid(), field::fzyx, FieldGhostLayers );
 
    // add boundary handling & initialize outer domain boundaries
    BlockDataID boundaryHandlingID = blocks->addBlockData( make_shared< MyBoundaryHandling >( blocks, flagFieldID, pdfFieldID, particleFieldID, accessor ), "boundary handling" );

apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadEvaluation.cpp

@@ -573,7 +573,7 @@ int main( int argc, char **argv )
 
       if(maxPenetrationDepth < overlapLimit) break;
 
-      // reset velocites to avoid too large ones
+      // reset velocities to avoid too large ones
 
       ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
             [](const size_t idx, ParticleAccessor_T& ac){
@@ -610,15 +610,15 @@ int main( int argc, char **argv )
    LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ) );
 
    // add PDF field
-   BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
+   BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (fzyx)", latticeModel,
                                                                          Vector3< real_t >( real_t(0) ), real_t(1),
-                                                                         uint_t(1), field::zyxf );
+                                                                         uint_t(1), field::fzyx );
 
    // add flag field
    BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field" );
 
    // add particle field
-   BlockDataID particleFieldID = field::addToStorage<lbm_mesapd_coupling::ParticleField_T>( blocks, "particle field", accessor->getInvalidUid(), field::zyxf, FieldGhostLayers );
+   BlockDataID particleFieldID = field::addToStorage<lbm_mesapd_coupling::ParticleField_T>( blocks, "particle field", accessor->getInvalidUid(), field::fzyx, FieldGhostLayers );
 
    // add boundary handling & initialize outer domain boundaries
    using BoundaryHandling_T = MyBoundaryHandling<ParticleAccessor_T>::Type;