diff --git a/apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp b/apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp new file mode 100644 index 0000000000000000000000000000000000000000..597b5ce7c3717c114cd64a0ffa255324d5a053b1 --- /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 d8234d517a84aec87c69e227dc875a3491f07cb7..19c2e838d94ffb69f5f2ecbfb93b037265c331cb 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 0000000000000000000000000000000000000000..ea8c5976585ea1c34a089828d6389e6f32ef1f8a --- /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)