refined SRT lattice model tutorial

apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp

@@ -28,6 +28,7 @@
 
 #include "geometry/all.h"
 
+#include "lbm/boundary/factories/DefaultBoundaryHandling.h"
 #include "lbm/communication/PdfFieldPackInfo.h"
 #include "lbm/field/AddToStorage.h"
 #include "lbm/field/PdfField.h"
@@ -37,57 +38,32 @@
 #include "timeloop/all.h"
 
 //    Codegen Includes
-#include "CumulantMRTLatticeModel.h"
-#include "EntropicMRTLatticeModel.h"
 #include "SRTLatticeModel.h"
-
-//    Generated Communication Pack Infos
-#include "CumulantMRTPackInfo.h"
-#include "EntropicMRTPackInfo.h"
 #include "SRTPackInfo.h"
 
-//    Generated NoSlip Boundaries
-#include "CumulantMRTNoSlip.h"
-#include "EntropicMRTNoSlip.h"
-#include "SRTNoSlip.h"
-
 namespace walberla
 {
 ///////////////////////
 /// Typedef Aliases ///
 ///////////////////////
 
-// Set a typedef alias for our generated lattice model and the boundary.
-// Changing the Method only requires changing these aliases.
-
+// Typedef alias for the lattice model
 typedef lbm::SRTLatticeModel LatticeModel_T;
-typedef lbm::SRTNoSlip NoSlip_T;
-
-// Entropic Model
-// typedef lbm::EntropicMRTLatticeModel LatticeModel_T;
-// typedef lbm::EntropicMRTNoSlip NoSlip_T;
-
-// Cumulant Model
-// typedef lbm::CumulantMRTLatticeModel LatticeModel_T;
-// typedef lbm::CumulantMRTNoSlip NoSlip_T;
 
 // Communication Pack Info
-typedef lbm::PdfFieldPackInfo< LatticeModel_T > PackInfo_T;
+typedef pystencils::SRTPackInfo PackInfo_T;
 
 // Also set aliases for the stencils involved...
 typedef LatticeModel_T::Stencil Stencil_T;
 typedef LatticeModel_T::CommunicationStencil CommunicationStencil_T;
 
-// ... and for the required field types. These include the PdfField for the simulation
+// ... and for the required field types.
 typedef lbm::PdfField< LatticeModel_T > PdfField_T;
 
-// and scalar- and vector-valued fields for output of macroscopic quantities.
-typedef GhostLayerField< real_t, LatticeModel_T::Stencil::D > VectorField_T;
-typedef GhostLayerField< real_t, 1 > ScalarField_T;
-
 // Also, for boundary handling, a flag data type and flag field are required.
 typedef walberla::uint8_t flag_t;
 typedef FlagField< flag_t > FlagField_T;
+typedef lbm::DefaultBoundaryHandlingFactory< LatticeModel_T, FlagField_T > BHFactory;
 
 /////////////////////////////////////////
 /// Shear Flow Initialization Functor ///
@@ -146,15 +122,8 @@ int main(int argc, char** argv)
    /// Field Setup ///
    ///////////////////
 
-   // TODO: Do we need a force field?
-   // BlockDataID forceFieldId = field::addToStorage<VectorField_T>( blocks, "Force", real_t( 0.0 ));
-
-   // Velocity output field
-   // We don't need that either, do we?
-   BlockDataID velFieldId = field::addToStorage< VectorField_T >(blocks, "Velocity", real_t(0.0));
-
-   LatticeModel_T latticeModel = LatticeModel_T(velFieldId, omega);
-   BlockDataID pdfFieldId      = lbm::addPdfFieldToStorage(blocks, "pdf field", latticeModel, field::zyxf);
+   LatticeModel_T latticeModel = LatticeModel_T(omega);
+   BlockDataID pdfFieldId      = lbm::addPdfFieldToStorage(blocks, "pdf field", latticeModel, field::fzyx);
    BlockDataID flagFieldId     = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
 
    ////////////////////////
@@ -174,12 +143,12 @@ int main(int argc, char** argv)
 
    auto boundariesConfig = walberlaEnv.config()->getOneBlock("Boundaries");
 
-   NoSlip_T noSlip(blocks, pdfFieldId);
+   BlockDataID boundaryHandlingId =
+      BHFactory::addBoundaryHandlingToStorage(blocks, "boundary handling", flagFieldId, pdfFieldId, fluidFlagUID,
+                                              Vector3< real_t >(), Vector3< real_t >(), real_c(0.0), real_c(0.0));
 
-   geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldId, boundariesConfig);
-   geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID);
-
-   noSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("NoSlip"), fluidFlagUID);
+   geometry::initBoundaryHandling< BHFactory::BoundaryHandling >(*blocks, boundaryHandlingId, boundariesConfig);
+   geometry::setNonBoundaryCellsToDomain< BHFactory::BoundaryHandling >(*blocks, boundaryHandlingId);
 
    /////////////////
    /// Time Loop ///
@@ -190,11 +159,11 @@ int main(int argc, char** argv)
    // Communication
    blockforest::communication::UniformBufferedScheme< CommunicationStencil_T > communication(blocks);
    communication.addPackInfo(make_shared< PackInfo_T >(pdfFieldId));
-   timeloop.add() << BeforeFunction(communication, "communication");
-
-   // Boundaries and LBM Sweep
-   timeloop.add() << Sweep(noSlip, "NoSlip Boundaries")
-                  << Sweep(LatticeModel_T::Sweep(pdfFieldId), "LB stream & collide");
+   
+   // Timeloop
+   timeloop.add() << BeforeFunction(communication, "communication")
+                  << Sweep(BHFactory::BoundaryHandling::getBlockSweep(boundaryHandlingId), "Boundary Handling");
+   timeloop.add() << Sweep(LatticeModel_T::Sweep(pdfFieldId), "LB stream & collide");
 
    // Stability Checker
    timeloop.addFuncAfterTimeStep(makeSharedFunctor(field::makeStabilityChecker< PdfField_T, FlagField_T >(

apps/tutorials/codegen/02_LBMLatticeModelGeneration.prm

@@ -35,7 +35,7 @@ StabilityChecker
 
 Boundaries 
 {   
-	Border { direction S,N; walldistance -1; flag NoSlip; }		
+	Border { direction S,N; walldistance -1; NoSlip {} }		
 }
 
 

apps/tutorials/codegen/02_LBMLatticeModelGeneration.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
+
+#   ========================
+#      General Parameters
+#   ========================
+
+STENCIL = 'D2Q9'
+OMEGA = sp.Symbol('omega')
+LAYOUT = 'fzyx'
+
+#   Optimization
+OPT = {'target': 'cpu', 'cse_global': True, 'field_layout': LAYOUT}
+
+#   ===========================
+#      SRT Method Definition
+#   ===========================
+
+srt_params = {'stencil': STENCIL,
+              'method': 'srt',
+              'relaxation_rate': OMEGA}
+
+srt_collision_rule = create_lb_collision_rule(optimization=OPT, **srt_params)
+srt_update_rule = create_lb_update_rule(collision_rule=srt_collision_rule, optimization=OPT)
+
+#   =====================
+#      Code Generation
+#   =====================
+
+with CodeGeneration() as ctx:
+    generate_lattice_model(ctx, "SRTLatticeModel", srt_collision_rule, field_layout=LAYOUT)
+    generate_pack_info_from_kernel(ctx, "SRTPackInfo", srt_update_rule)

apps/tutorials/codegen/CMakeLists.txt

@@ -12,22 +12,15 @@ if( WALBERLA_BUILD_WITH_CODEGEN )
                               DEPENDS blockforest core field stencil timeloop vtk pde CodegenHeatEquationKernel )
 
     #   Tutorial 2: lbmpy Lattice Model Generation
-    walberla_generate_target_from_python( NAME LBMLatticeModelGenerationPython
-        FILE LBMLatticeModelGeneration.py
-        OUT_FILES   SRTLatticeModel.cpp SRTLatticeModel.h 
-                    SRTPackInfo.cpp SRTPackInfo.h
-                    SRTNoSlip.cpp SRTNoSlip.h
-                    CumulantMRTLatticeModel.cpp CumulantMRTLatticeModel.h 
-                    CumulantMRTPackInfo.cpp CumulantMRTPackInfo.h
-                    CumulantMRTNoSlip.cpp CumulantMRTNoSlip.h
-                    EntropicMRTLatticeModel.cpp EntropicMRTLatticeModel.h
-                    EntropicMRTPackInfo.cpp EntropicMRTPackInfo.h
-                    EntropicMRTNoSlip.cpp EntropicMRTNoSlip.h )
-
     waLBerla_link_files_to_builddir( *.prm )
+    
+    walberla_generate_target_from_python( NAME 02_LBMLatticeModelGenerationPython
+        FILE 02_LBMLatticeModelGeneration.py
+        OUT_FILES   SRTLatticeModel.cpp SRTLatticeModel.h 
+                    SRTPackInfo.cpp SRTPackInfo.h )
 
     walberla_add_executable ( NAME 02_LBMLatticeModelGenerationApp 
                     FILES 02_LBMLatticeModelGeneration.cpp
-                    DEPENDS blockforest core domain_decomposition field geometry timeloop lbm stencil vtk LBMLatticeModelGenerationPython )
+                    DEPENDS blockforest core domain_decomposition field geometry timeloop lbm stencil vtk 02_LBMLatticeModelGenerationPython )
 
 endif()
\ No newline at end of file

apps/tutorials/codegen/LBMLatticeModelGeneration.py

-import sympy as sp
-import pystencils as ps
-
-from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_update_rule
-from lbmpy.boundaries import NoSlip
-
-from pystencils_walberla import CodeGeneration, generate_pack_info_from_kernel
-from lbmpy_walberla import generate_lattice_model, generate_boundary
-
-#   =====================
-#   Common Parameters
-#   =====================
-
-STENCIL = 'D2Q9'
-OMEGA = sp.Symbol('omega')
-LAYOUT = 'zyxf'
-
-#   Optimization
-OPT = {'target': 'cpu', 'cse_global': True, 'field_layout': LAYOUT}
-
-#   Velocity Output Field
-vel_field = ps.fields("velocity(2): [2D]", layout=LAYOUT)
-OUTPUT = {'velocity': vel_field}
-
-#   ==================================================================================================
-#   Method Definitions, including collision rules, boundary handling and communication pack infos.
-#   ==================================================================================================
-
-
-#   SRT Method
-def build_srt_method(ctx):
-    srt_params = {'stencil': STENCIL,
-                  'method': 'srt',
-                  'relaxation_rate': OMEGA}
-
-    srt_collision_rule = create_lb_collision_rule(optimization=OPT, output=OUTPUT, **srt_params)
-    generate_lattice_model(ctx, "SRTLatticeModel", srt_collision_rule, field_layout=LAYOUT)
-
-    srt_update_rule = create_lb_update_rule(collision_rule=srt_collision_rule, optimization=OPT)
-    generate_pack_info_from_kernel(ctx, "SRTPackInfo", srt_update_rule)
-
-    generate_boundary(ctx, "SRTNoSlip", NoSlip(), srt_collision_rule.method)
-
-
-#   Cumulant MRT Method
-def build_cumulant_method(ctx):
-    mrt_cumulant_params = {'stencil': STENCIL,
-                           'method': 'mrt_raw',
-                           'relaxation_rates': [0, 0, 0, OMEGA, OMEGA, OMEGA, 1, 1, 1],
-                           'cumulant': True,
-                           'compressible': True}   # Incompressible cumulants not yet supported!
-
-    mrt_cumulant_collision_rule = create_lb_collision_rule(optimization=OPT, output=OUTPUT, **mrt_cumulant_params)
-    generate_lattice_model(ctx, "CumulantMRTLatticeModel", mrt_cumulant_collision_rule, field_layout=LAYOUT)
-
-    mrt_cumulant_update_rule = create_lb_update_rule(collision_rule=mrt_cumulant_collision_rule, optimization=OPT)
-    generate_pack_info_from_kernel(ctx, "CumulantMRTPackInfo", mrt_cumulant_update_rule)
-
-    generate_boundary(ctx, "CumulantMRTNoSlip", NoSlip(), mrt_cumulant_collision_rule.method)
-
-
-#   Orthogonal MRT Method with entropy constraint
-def build_entropic_method(ctx):
-    omega_f = sp.Symbol('omega_f')
-
-    mrt_entropic_params = {'stencil': STENCIL,
-                           'method': 'mrt',
-                           'relaxation_rates': [OMEGA, OMEGA, omega_f, omega_f],
-                           'entropic': True,
-                           'compressible': True}    # Entropic models only implemented with pdfs centered around 1
-
-    mrt_entropic_collision_rule = create_lb_collision_rule(optimization=OPT, output=OUTPUT, **mrt_entropic_params)
-    generate_lattice_model(ctx, "EntropicMRTLatticeModel", mrt_entropic_collision_rule, field_layout=LAYOUT)
-
-    mrt_entropic_update_rule = create_lb_update_rule(collision_rule=mrt_entropic_collision_rule, optimization=OPT)
-    generate_pack_info_from_kernel(ctx, "EntropicMRTPackInfo", mrt_entropic_update_rule)
-
-    generate_boundary(ctx, "EntropicMRTNoSlip", NoSlip(), mrt_entropic_collision_rule.method)
-
-#   ================================
-#   Main Function
-#   ================================
-
-
-if __name__ == "__main__":
-    #   All code generation must happen within one context, and all files specified in
-    #   CMakeLists.txt must be generated therein.
-    with CodeGeneration() as ctx:
-        build_srt_method(ctx)
-        build_cumulant_method(ctx)
-        build_entropic_method(ctx)