Add minimal example of channel flow for comparing walberla and lbmpy SRT

apps/CMakeLists.txt

@@ -32,4 +32,6 @@ endif()
 # Python module
 if ( WALBERLA_BUILD_WITH_PYTHON )
     add_subdirectory( pythonmodule )
-endif()
\ No newline at end of file
+endif()
+
+add_subdirectory( channelFlow )
\ No newline at end of file

apps/channelFlow/CMakeLists.txt

+waLBerla_link_files_to_builddir( *.prm )
+
+waLBerla_add_executable ( NAME walberlaSRT
+                          FILES walberlaSRT.cpp
+                          DEPENDS blockforest core field lbm geometry timeloop gui )
+
+if( WALBERLA_BUILD_WITH_CODEGEN )
+    walberla_generate_target_from_python( NAME srtLatticeModel
+                                          FILE latticeModel.py
+                                          OUT_FILES   SRTLatticeModel.cpp SRTLatticeModel.h )
+
+    waLBerla_add_executable ( NAME lbmpySRT
+                              FILES lbmpySRT.cpp
+                              DEPENDS blockforest core field lbm geometry timeloop gui srtLatticeModel)
+endif()
\ No newline at end of file

apps/channelFlow/channelFlow.prm

+
+Parameters 
+{
+	omega           1.8;
+	initialVelocity < 0.001, 0, 0 >;
+	timesteps       10;
+}
+
+DomainSetup
+{
+   blocks        <  1,    1, 1 >;
+   cellsPerBlock <  5, 5, 1 >;
+   periodic      <  1,    0, 1 >;
+}
+
+Boundaries 
+{
+	Border { direction S,N;  walldistance -1;  NoSlip    {} }
+}
+
+VTK 
+{
+   fluid_field
+   {
+      writeFrequency 1;
+      
+      inclusion_filters {
+         DomainFilter;
+      }
+      
+      writers {
+         Velocity;
+         Density;
+      }
+   }
+   
+   flag_field
+   {
+      writeFrequency 10000000; // write only once
+      
+      writers {
+         FlagField;
+      }
+   }
+}

apps/channelFlow/latticeModel.py

+import sympy as sp
+
+from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_update_rule
+
+from pystencils_walberla import CodeGeneration, generate_pack_info_from_kernel
+from lbmpy_walberla import generate_lattice_model
+
+stencil = 'D2Q9'
+omega = sp.Symbol('omega')
+layout = 'fzyx'
+
+optimizations = {'cse_global': True, 'field_layout': layout}
+
+srt_params = {'stencil': stencil,
+              'method': 'srt',
+              'relaxation_rate': omega,
+              'compressible': True,
+              'equilibrium_order': 2}
+
+srt_collision_rule = create_lb_collision_rule(optimization=optimizations, **srt_params)
+srt_update_rule = create_lb_update_rule(collision_rule=srt_collision_rule, optimization=optimizations)
+
+with CodeGeneration() as ctx:
+    generate_lattice_model(ctx, "SRTLatticeModel", srt_collision_rule, field_layout=layout)

apps/channelFlow/lbmpySRT.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file lbmpySRT.cpp
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/all.h"
+
+#include "core/all.h"
+
+#include "domain_decomposition/all.h"
+
+#include "field/all.h"
+
+#include "geometry/all.h"
+
+#include "lbm/all.h"
+
+#include "timeloop/all.h"
+
+#include "SRTLatticeModel.h"
+
+namespace walberla
+{
+using LatticeModel_T         = lbm::SRTLatticeModel;
+using Stencil_T              = LatticeModel_T::Stencil;
+using CommunicationStencil_T = LatticeModel_T::CommunicationStencil;
+
+using PdfField_T = lbm::PdfField< LatticeModel_T >;
+
+using flag_t      = walberla::uint8_t;
+using FlagField_T = FlagField< flag_t >;
+
+int main(int argc, char** argv)
+{
+   walberla::Environment walberlaEnv(argc, argv);
+
+   auto blocks = blockforest::createUniformBlockGridFromConfig(walberlaEnv.config());
+
+   // read parameters
+   auto parameters = walberlaEnv.config()->getOneBlock("Parameters");
+
+   const real_t omega = parameters.getParameter< real_t >("omega", real_c(1.4));
+   const Vector3< real_t > initialVelocity =
+      parameters.getParameter< Vector3< real_t > >("initialVelocity", Vector3< real_t >());
+   const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
+
+   // create fields
+   LatticeModel_T latticeModel = LatticeModel_T(omega);
+   BlockDataID pdfFieldId =
+      lbm::addPdfFieldToStorage(blocks, "pdf field", latticeModel, initialVelocity, real_t(1), field::fzyx);
+   BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+
+   // create and initialize boundary handling
+   const FlagUID fluidFlagUID("Fluid");
+
+   auto boundariesConfig = walberlaEnv.config()->getOneBlock("Boundaries");
+
+   using BHFactory = lbm::DefaultBoundaryHandlingFactory< LatticeModel_T, FlagField_T >;
+
+   BlockDataID boundaryHandlingId = BHFactory::addBoundaryHandlingToStorage(
+      blocks, "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.0)));
+
+   geometry::initBoundaryHandling< BHFactory::BoundaryHandling >(*blocks, boundaryHandlingId, boundariesConfig);
+   geometry::setNonBoundaryCellsToDomain< BHFactory::BoundaryHandling >(*blocks, boundaryHandlingId);
+
+   // create time loop
+   SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
+
+   // create communication for PdfField
+   blockforest::communication::UniformBufferedScheme< CommunicationStencil_T > communication(blocks);
+   communication.addPackInfo(make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >(pdfFieldId));
+
+   // add LBM sweep and communication to time loop
+   timeloop.add() << BeforeFunction(communication, "communication")
+                  << Sweep(BHFactory::BoundaryHandling::getBlockSweep(boundaryHandlingId), "boundary handling");
+
+   auto generatedLbmSweep = make_shared< typename LatticeModel_T::Sweep >(pdfFieldId);
+   timeloop.add() << Sweep(lbm::makeStreamSweep(generatedLbmSweep), "Stream")
+                  << Sweep(lbm::makeCollideSweep(generatedLbmSweep), "Collide");
+   // timeloop.add() << Sweep(*generatedLbmSweep, "LB stream & collide");
+
+   // add VTK output to time loop
+   lbm::VTKOutput< LatticeModel_T, FlagField_T >::addToTimeloop(timeloop, blocks, walberlaEnv.config(), pdfFieldId,
+                                                                flagFieldId, fluidFlagUID);
+
+   timeloop.run();
+
+   return EXIT_SUCCESS;
+}
+} // namespace walberla
+
+int main(int argc, char** argv) { walberla::main(argc, argv); }
\ No newline at end of file

apps/channelFlow/walberlaSRT.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file walberlaSRT.cpp
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/all.h"
+
+#include "core/all.h"
+
+#include "domain_decomposition/all.h"
+
+#include "field/all.h"
+
+#include "geometry/all.h"
+
+#include "lbm/all.h"
+
+#include "timeloop/all.h"
+
+namespace walberla
+{
+using LatticeModel_T         = lbm::D2Q9< lbm::collision_model::SRT, true, lbm::force_model::None, 2 >;
+using Stencil_T              = LatticeModel_T::Stencil;
+using CommunicationStencil_T = LatticeModel_T::CommunicationStencil;
+
+using PdfField_T = lbm::PdfField< LatticeModel_T >;
+
+using flag_t      = walberla::uint8_t;
+using FlagField_T = FlagField< flag_t >;
+
+int main(int argc, char** argv)
+{
+   walberla::Environment walberlaEnv(argc, argv);
+
+   auto blocks = blockforest::createUniformBlockGridFromConfig(walberlaEnv.config());
+
+   // read parameters
+   auto parameters = walberlaEnv.config()->getOneBlock("Parameters");
+
+   const real_t omega = parameters.getParameter< real_t >("omega", real_c(1.4));
+   const Vector3< real_t > initialVelocity =
+      parameters.getParameter< Vector3< real_t > >("initialVelocity", Vector3< real_t >());
+   const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
+
+   // create fields
+   LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::SRT(omega));
+   BlockDataID pdfFieldId =
+      lbm::addPdfFieldToStorage(blocks, "pdf field", latticeModel, initialVelocity, real_t(1), field::fzyx);
+   BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+
+   // create and initialize boundary handling
+   const FlagUID fluidFlagUID("Fluid");
+
+   auto boundariesConfig = walberlaEnv.config()->getOneBlock("Boundaries");
+
+   using BHFactory = lbm::DefaultBoundaryHandlingFactory< LatticeModel_T, FlagField_T >;
+
+   BlockDataID boundaryHandlingId = BHFactory::addBoundaryHandlingToStorage(
+      blocks, "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.0)));
+
+   geometry::initBoundaryHandling< BHFactory::BoundaryHandling >(*blocks, boundaryHandlingId, boundariesConfig);
+   geometry::setNonBoundaryCellsToDomain< BHFactory::BoundaryHandling >(*blocks, boundaryHandlingId);
+
+   // create time loop
+   SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
+
+   // create communication for PdfField
+   blockforest::communication::UniformBufferedScheme< CommunicationStencil_T > communication(blocks);
+   communication.addPackInfo(make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >(pdfFieldId));
+
+   // add LBM sweep and communication to time loop
+   timeloop.add() << BeforeFunction(communication, "communication")
+                  << Sweep(BHFactory::BoundaryHandling::getBlockSweep(boundaryHandlingId), "boundary handling");
+   timeloop.add() << Sweep(
+      makeSharedSweep(lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >(pdfFieldId, flagFieldId, fluidFlagUID)),
+      "LB stream & collide");
+
+   // add VTK output to time loop
+   lbm::VTKOutput< LatticeModel_T, FlagField_T >::addToTimeloop(timeloop, blocks, walberlaEnv.config(), pdfFieldId,
+                                                                flagFieldId, fluidFlagUID);
+
+   timeloop.run();
+
+   return EXIT_SUCCESS;
+}
+} // namespace walberla
+
+int main(int argc, char** argv) { walberla::main(argc, argv); }
\ No newline at end of file