From 68c49ac95b17b69898813b21e754c42154f119fe Mon Sep 17 00:00:00 2001
From: Frederik Hennig <frederik.hennig@fau.de>
Date: Wed, 26 Aug 2020 11:43:20 +0200
Subject: [PATCH] lbmpy codegen script, stub application source, and
CMakeLists.
---
.../codegen/02_LBMLatticeModelGeneration.cpp | 132 ++++++++++++++++++
apps/tutorials/codegen/CMakeLists.txt | 20 +++
.../codegen/LBMLatticeModelGeneration.py | 91 ++++++++++++
3 files changed, 243 insertions(+)
create mode 100644 apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp
create mode 100644 apps/tutorials/codegen/LBMLatticeModelGeneration.py
diff --git a/apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp b/apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp
new file mode 100644
index 000000000..597b5ce7c
--- /dev/null
+++ b/apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp
@@ -0,0 +1,132 @@
+//======================================================================================================================
+//
+// 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 02_LBMLatticeModelGeneration.cpp
+//! \author Frederik Hennig <frederik.hennig@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/all.h"
+#include "core/all.h"
+#include "domain_decomposition/all.h"
+#include "field/all.h"
+#include "geometry/all.h"
+#include "timeloop/all.h"
+
+#include "lbm/field/PdfField.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/vtk/VTKOutput.h"
+
+// Codegen Includes
+#include "SRTLatticeModel.h"
+#include "CumulantMRTLatticeModel.h"
+#include "EntropicMRTLatticeModel.h"
+
+namespace walberla{
+
+// Set a typedef alias for our generated lattice model, for convenience.
+typedef lbm::SRTLatticeModel LatticeModel_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
+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;
+
+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 ) );
+
+ const double remainingTimeLoggerFrequency = parameters.getParameter< double >( "remainingTimeLoggerFrequency", 3.0 ); // in seconds
+
+ // create fields
+
+ // TODO: Do we need a force field?
+ //BlockDataID forceFieldId = field::addToStorage<VectorField_T>( blocks, "Force", real_t( 0.0 ));
+
+ // Velocity output field
+ 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, initialVelocity, real_t(1) );
+ BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
+
+ // create and initialize boundary handling
+ const FlagUID fluidFlagUID( "Fluid" );
+
+
+ auto boundariesConfig = walberlaEnv.config()->getOneBlock( "Boundaries" );
+
+ // TODO: Replace by my NoSlip
+ // lbm::LbCodeGenerationExample_NoSlip noSlip(blocks, pdfFieldId);
+
+ geometry::initBoundaryHandling<FlagField_T>(*blocks, flagFieldId, boundariesConfig);
+ geometry::setNonBoundaryCellsToDomain<FlagField_T>(*blocks, flagFieldId, fluidFlagUID);
+
+ // TODO
+ // noSlip.fillFromFlagField<FlagField_T>( blocks, flagFieldId, FlagUID("NoSlip"), fluidFlagUID );
+
+ // 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( walls, "Boundary Walls" );
+ timeloop.add() << Sweep( LatticeModel_T::Sweep( pdfFieldId ), "LB stream & collide" );
+
+ // LBM stability check
+ timeloop.addFuncAfterTimeStep( makeSharedFunctor( field::makeStabilityChecker< PdfField_T, FlagField_T >( walberlaEnv.config(), blocks, pdfFieldId,
+ flagFieldId, fluidFlagUID ) ),
+ "LBM stability check" );
+
+ // log remaining time
+ timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger( timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency ), "remaining time logger" );
+
+ // 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;
+}
+
+}
+
+int main(int argc, char ** argv){
+ return walberla::main(argc, argv);
+}
\ No newline at end of file
diff --git a/apps/tutorials/codegen/CMakeLists.txt b/apps/tutorials/codegen/CMakeLists.txt
index d8234d517..19c2e838d 100644
--- a/apps/tutorials/codegen/CMakeLists.txt
+++ b/apps/tutorials/codegen/CMakeLists.txt
@@ -1,6 +1,8 @@
# Code Generation Tutorials
if( WALBERLA_BUILD_WITH_CODEGEN )
+
+ # Tutorial 1: Heat Equation with pystencils
walberla_generate_target_from_python( NAME CodegenHeatEquationKernel
FILE HeatEquationKernel.py
OUT_FILES HeatEquationKernel.cpp HeatEquationKernel.h )
@@ -8,4 +10,22 @@ if( WALBERLA_BUILD_WITH_CODEGEN )
walberla_add_executable ( NAME 01_CodegenHeatEquation
FILES 01_CodegenHeatEquation.cpp
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_add_executable ( NAME 02_LBMLatticeModelGenerationApp
+ FILES 02_LBMLatticeModelGeneration.cpp
+ DEPENDS blockforest core field stencil timeloop vtk pde LBMLatticeModelGenerationPython )
+
endif()
\ No newline at end of file
diff --git a/apps/tutorials/codegen/LBMLatticeModelGeneration.py b/apps/tutorials/codegen/LBMLatticeModelGeneration.py
new file mode 100644
index 000000000..ea8c59765
--- /dev/null
+++ b/apps/tutorials/codegen/LBMLatticeModelGeneration.py
@@ -0,0 +1,91 @@
+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}
+
+# Velocity Output Field
+vel_field = ps.fields("velocity(3): [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)
+ 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)
+ 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)
+ 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)
--
GitLab