Added ChannelFlow

apps/showcases/CMakeLists.txt

@@ -7,7 +7,10 @@ add_subdirectory( Mixer )
 add_subdirectory( ParticlePacking )
 add_subdirectory( PegIntoSphereBed )
 if ( WALBERLA_BUILD_WITH_CODEGEN AND WALBERLA_BUILD_WITH_PYTHON )
-add_subdirectory( PhaseFieldAllenCahn )
+   add_subdirectory( ChannelFlow )
+   if (WALBERLA_BUILD_WITH_PYTHON)
+      add_subdirectory( PhaseFieldAllenCahn )
+   endif()
 endif()
 if ( WALBERLA_BUILD_WITH_CODEGEN AND NOT WALBERLA_BUILD_WITH_OPENMP)
    add_subdirectory( PorousMedia )

apps/showcases/ChannelFlow/CMakeLists.txt

+waLBerla_link_files_to_builddir( "*.prm" )
+
+waLBerla_generate_target_from_python(NAME ChannelFlowGenerated
+        FILE ChannelFlow.py
+        OUT_FILES ChannelFlow_Sweep.cpp ChannelFlow_Sweep.h
+        ChannelFlow_MacroSetter.cpp ChannelFlow_MacroSetter.h
+        ChannelFlow_NoSlip.cpp ChannelFlow_NoSlip.h
+        ChannelFlow_PackInfo.cpp ChannelFlow_PackInfo.h
+        ChannelFlow.h)
+waLBerla_add_executable(NAME ChannelFlow
+        FILES ChannelFlow.cpp
+        DEPENDS blockforest core field lbm geometry timeloop ChannelFlowGenerated)

apps/showcases/ChannelFlow/ChannelFlow.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 ChannelFlow.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "core/Environment.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "field/AddToStorage.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "geometry/InitBoundaryHandling.h"
+#include "timeloop/SweepTimeloop.h"
+
+// Generated Files
+#include "ChannelFlow.h"
+
+using namespace walberla;
+
+using PackInfo_T  = lbm::ChannelFlow_PackInfo;
+using flag_t      = walberla::uint8_t;
+using FlagField_T = FlagField< flag_t >;
+
+auto pdfFieldAdder = [](IBlock* const block, StructuredBlockStorage* const storage) {
+   return new PdfField_T(storage->getNumberOfXCells(*block), storage->getNumberOfYCells(*block),
+                         storage->getNumberOfZCells(*block), uint_t(1), field::fzyx,
+                         make_shared< field::AllocateAligned< real_t, 64 > >());
+};
+
+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");
+   const uint_t timesteps = parameters.getParameter< uint_t >("timesteps");
+   real_t ForceX = parameters.getParameter< real_t >("ForceX");
+
+   WALBERLA_LOG_INFO_ON_ROOT("Using a relaxation rate of " << omega)
+   WALBERLA_LOG_INFO_ON_ROOT("Set force " << ForceX << " in x direction.")
+
+   const double remainingTimeLoggerFrequency =
+      parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
+
+   // create fields
+   BlockDataID pdfFieldID     = blocks->addStructuredBlockData< PdfField_T >(pdfFieldAdder, "PDFs");
+   BlockDataID velFieldID     = field::addToStorage< VelocityField_T >(blocks, "velocity", real_t(0.0), field::fzyx);
+   BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(1.0), field::fzyx);
+
+   BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+
+   pystencils::ChannelFlow_MacroSetter setterSweep(densityFieldID, pdfFieldID, velFieldID, ForceX);
+   for (auto& block : *blocks)
+      setterSweep(&block);
+
+   // create and initialize boundary handling
+   const FlagUID fluidFlagUID("Fluid");
+
+   auto boundariesConfig = walberlaEnv.config()->getOneBlock("Boundaries");
+
+   lbm::ChannelFlow_NoSlip noSlip(blocks, pdfFieldID);
+
+   geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldId, boundariesConfig);
+   geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID);
+
+   noSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("NoSlip"), fluidFlagUID);
+
+   // create time loop
+   SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
+
+   // create communication for PdfField
+   blockforest::communication::UniformBufferedScheme< Stencil_T > communication(blocks);
+   communication.addPackInfo(make_shared< PackInfo_T >(pdfFieldID));
+
+   pystencils::ChannelFlow_Sweep UpdateSweep(densityFieldID, pdfFieldID, velFieldID, ForceX, omega);
+
+   // add LBM sweep and communication to time loop
+   timeloop.add() << Sweep(noSlip, "noSlip boundary") << AfterFunction(communication, "communication");
+   timeloop.add() << Sweep(UpdateSweep, "LB stream & collide");
+
+   // log remaining time
+   timeloop.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
+                                 "remaining time logger");
+
+   // VTK Writer
+   uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0);
+   if (vtkWriteFrequency > 0)
+   {
+      auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "ChannelFlow_VTK", vtkWriteFrequency, 0, false,
+                                                      "vtk_out", "simulation_step", false, true, true, false, 0);
+
+      auto velWriter     = make_shared< field::VTKWriter< VelocityField_T > >(velFieldID, "velocity");
+      auto densityWriter = make_shared< field::VTKWriter< ScalarField_T > >(densityFieldID, "density");
+
+      vtkOutput->addCellDataWriter(velWriter);
+      vtkOutput->addCellDataWriter(densityWriter);
+
+      timeloop.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
+   }
+
+   WALBERLA_LOG_INFO_ON_ROOT("Starting Simulation with " << timesteps << " timesteps")
+   timeloop.run();
+   WALBERLA_LOG_INFO_ON_ROOT("Finished Simulation")
+
+   return EXIT_SUCCESS;
+}

apps/showcases/ChannelFlow/ChannelFlow.prm

+
+Parameters 
+{
+	omega           1.99808;
+	timesteps       5000001;
+	ForceX          1.5e-7;
+	vtkWriteFrequency 100000;
+	remainingTimeLoggerFrequency 30; // in seconds
+}
+
+DomainSetup
+{
+   blocks        <  1,    1, 1 >;
+   cellsPerBlock <  48, 32, 16 >;
+   periodic      <  1,    0, 1 >;
+}
+
+Boundaries 
+{
+    Border { direction N;    walldistance -1;  flag NoSlip; }
+    Border { direction S;    walldistance -1;  flag NoSlip; }
+}

apps/showcases/ChannelFlow/ChannelFlow.py

+import sympy as sp
+
+from pystencils.field import fields
+from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
+from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Stencil, Method
+from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
+from lbmpy.boundaries import NoSlip
+from pystencils_walberla import CodeGeneration, generate_sweep, generate_info_header
+from lbmpy_walberla import generate_boundary, generate_lb_pack_info
+
+
+stencil = LBStencil(Stencil.D3Q19)
+
+pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): double[{stencil.D}D]", layout='fzyx')
+velocity_field, density_field = fields(f"velocity({stencil.D}), density(1) : double[{stencil.D}D]", layout='fzyx')
+omega = sp.Symbol("omega")
+fx = sp.symbols("fx")
+
+output = {
+    'density': density_field,
+    'velocity': velocity_field
+}
+
+lbm_config = LBMConfig(method=Method.SRT, stencil=stencil, relaxation_rate=omega, force=(fx, 0, 0),
+                       compressible=False,
+                       field_name='pdfs', output=output)
+
+lbm_opt = LBMOptimisation(symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp,
+                          cse_global=False, cse_pdfs=False)
+
+method = create_lb_method(lbm_config=lbm_config)
+
+# getter & setter
+setter_assignments = macroscopic_values_setter(method, velocity=velocity_field.center_vector,
+                                               pdfs=pdfs, density=density_field)
+
+update_rule = create_lb_update_rule(lb_method=method, lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+
+stencil_typedefs = {'Stencil_T': stencil}
+field_typedefs = {'PdfField_T': pdfs,
+                  'VelocityField_T': velocity_field,
+                  'ScalarField_T': density_field}
+
+with CodeGeneration() as ctx:
+    # sweeps
+    generate_sweep(ctx, 'ChannelFlow_Sweep', update_rule, field_swaps=[(pdfs, pdfs_tmp)])
+    generate_sweep(ctx, 'ChannelFlow_MacroSetter', setter_assignments)
+
+    generate_boundary(ctx, 'ChannelFlow_NoSlip', NoSlip(), method)
+
+    # communication
+    generate_lb_pack_info(ctx, 'ChannelFlow_PackInfo', stencil, pdfs)
+
+    # Info header containing correct template definitions for stencil and field
+    generate_info_header(ctx, "ChannelFlow.h",
+                         stencil_typedefs=stencil_typedefs, field_typedefs=field_typedefs)