Merge branch 'Phasefield_Allen_Cahn' into 'master'

Phasefield allen cahn

See merge request walberla/walberla!302

AUTHORS.txt

 List of contributors
 ====================
 
+Cameron Stewart
 Christian Feichtinger
 Christian Godenschwager
 Christoph Rettinger
@@ -12,14 +13,20 @@ Dominik Bartuschat
 Ehsan Fattahi
 Felix Winterhalter
 Florian Schornbaum
+Helen Schottenhamml
 Jan Götz
+Jan Hönig
+João Victor Tozatti Risso
 Johannes Habich
 Klaus Iglberger
 Kristina Pickl
 Lorenz Hufnagel
+Lukas Werner
+Markus Holzer
 Martin Bauer
 Matthias Markl
 Michael Kuron
+Nils Kohl
 Paulo Carvalho
 Regina Ammer
 Sagar Dolas
@@ -27,7 +34,9 @@ Sebastian Eibl
 Silke Bergler
 Simon Bogner
 Stefan Donath
+Stephan Seitz
 Sunil Kontham
+Tobias Leemann
 Tobias Preclik
 Tobias Scharpff
 Tobias Schruff

apps/benchmarks/CMakeLists.txt

@@ -14,9 +14,11 @@ add_subdirectory( PoiseuilleChannel )
 add_subdirectory( ProbeVsExtraMessage )
 add_subdirectory( SchaeferTurek )
 add_subdirectory( UniformGrid )
-if ( WALBERLA_BUILD_WITH_CODEGEN )
+if ( WALBERLA_BUILD_WITH_CODEGEN AND NOT WALBERLA_BUILD_WITH_CUDA )
 add_subdirectory( UniformGridGenerated )
+add_subdirectory( PhaseFieldAllenCahn )
 endif()
 if ( WALBERLA_BUILD_WITH_CUDA )
 add_subdirectory( UniformGridGPU )
-endif()
+add_subdirectory( PhaseFieldAllenCahn )
\ No newline at end of file
+endif()

apps/benchmarks/PhaseFieldAllenCahn/CMakeLists.txt

+waLBerla_link_files_to_builddir(*.prm)
+waLBerla_link_files_to_builddir(*.py)
+
+if (WALBERLA_BUILD_WITH_CUDA)
+    waLBerla_generate_target_from_python(NAME BenchmarkPhaseFieldCodeGenGPU
+            FILE multiphase_codegen.py
+            OUT_FILES initialize_phase_field_distributions.cu initialize_phase_field_distributions.h
+            initialize_velocity_based_distributions.cu initialize_velocity_based_distributions.h
+            phase_field_LB_step.cu phase_field_LB_step.h
+            hydro_LB_step.cu hydro_LB_step.h
+            PackInfo_phase_field_distributions.cu PackInfo_phase_field_distributions.h
+            PackInfo_phase_field.cu PackInfo_phase_field.h
+            PackInfo_velocity_based_distributions.cu PackInfo_velocity_based_distributions.h
+            GenDefines.h)
+
+    waLBerla_add_executable(NAME benchmark_multiphase
+            FILES benchmark_multiphase.cpp InitializerFunctions.cpp multiphase_codegen.py
+            DEPENDS blockforest core cuda field postprocessing lbm geometry timeloop gui BenchmarkPhaseFieldCodeGenGPU)
+else ()
+    waLBerla_generate_target_from_python(NAME BenchmarkPhaseFieldCodeGenCPU
+            FILE multiphase_codegen.py
+            OUT_FILES initialize_phase_field_distributions.cpp initialize_phase_field_distributions.h
+            initialize_velocity_based_distributions.cpp initialize_velocity_based_distributions.h
+            phase_field_LB_step.cpp phase_field_LB_step.h
+            hydro_LB_step.cpp hydro_LB_step.h
+            PackInfo_phase_field_distributions.cpp PackInfo_phase_field_distributions.h
+            PackInfo_phase_field.cpp PackInfo_phase_field.h
+            PackInfo_velocity_based_distributions.cpp PackInfo_velocity_based_distributions.h
+            GenDefines.h)
+
+    waLBerla_add_executable(NAME benchmark_multiphase
+            FILES benchmark_multiphase.cpp InitializerFunctions.cpp multiphase_codegen.py
+            DEPENDS blockforest core field postprocessing lbm geometry timeloop gui BenchmarkPhaseFieldCodeGenCPU)
+endif (WALBERLA_BUILD_WITH_CUDA)
+

apps/benchmarks/PhaseFieldAllenCahn/InitializerFunctions.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 InitializerFunctions.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+#include "core/math/Constants.h"
+
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "python_coupling/DictWrapper.h"
+
+namespace walberla
+{
+using PhaseField_T = GhostLayerField< real_t, 1 >;
+
+void initPhaseField_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
+                           const real_t R                         = 10,
+                           const Vector3< real_t > bubbleMidPoint = Vector3< real_t >(0.0, 0.0, 0.0),
+                           const real_t W                         = 5)
+{
+   for (auto& block : *blocks)
+   {
+      auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
+      // clang-format off
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell;
+         blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+         real_t Ri = sqrt((globalCell[0] - bubbleMidPoint[0]) * (globalCell[0] - bubbleMidPoint[0]) +
+                          (globalCell[1] - bubbleMidPoint[1]) * (globalCell[1] - bubbleMidPoint[1]) +
+                          (globalCell[2] - bubbleMidPoint[2]) * (globalCell[2] - bubbleMidPoint[2]));
+         phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (Ri - R) / W);
+      )
+      // clang-format on
+   }
+}
+
+void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
+                        const real_t W = 5)
+{
+   auto X              = blocks->getDomainCellBB().xMax();
+   auto halfY          = (blocks->getDomainCellBB().yMax()) / 2.0;
+   double perturbation = 0.05;
+
+   for (auto& block : *blocks)
+   {
+      auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
+      // clang-format off
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell;
+         blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+         real_t tmp = perturbation * X * (cos((2.0 * math::pi * globalCell[0]) / X) + cos((2.0 * math::pi * globalCell[2]) / X));
+         phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) - tmp) / (W / 2.0));
+      )
+      // clang-format on
+   }
+}
+} // namespace walberla

apps/benchmarks/PhaseFieldAllenCahn/InitializerFunctions.h

+//======================================================================================================================
+//
+//  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 InitializerFunctions.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+#include "core/math/Constants.h"
+
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "python_coupling/DictWrapper.h"
+#pragma once
+
+namespace walberla
+{
+void initPhaseField_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
+                           Vector3< real_t > bubbleMidPoint, real_t W = 5);
+
+void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W = 5);
+
+} // namespace walberla

apps/benchmarks/PhaseFieldAllenCahn/benchmark.py

+import os
+import waLBerla as wlb
+import pandas as pd
+
+from waLBerla.tools.sqlitedb import sequenceValuesToScalars
+
+import sys
+
+
+class Scenario:
+    def __init__(self, timeStepStrategy, overlappingWidth, cuda_block_size):
+        # output frequencies
+        self.vtkWriteFrequency = 0
+
+        # simulation parameters
+        self.timesteps = 101
+        edge_size = 32
+        self.cells = (edge_size, edge_size, edge_size)
+        self.blocks = (1, 1, 1)
+        self.periodic = (1, 1, 1)
+        self.size = (self.cells[0] * self.blocks[0],
+                     self.cells[1] * self.blocks[1],
+                     self.cells[2] * self.blocks[2])
+
+        self.timeStepStrategy = timeStepStrategy
+        self.overlappingWidth = overlappingWidth
+        self.cuda_block_size = cuda_block_size
+        self.warmupSteps = 10
+
+        # bubble parameters
+        self.bubbleRadius = min(self.size) // 4
+        self.bubbleMidPoint = (self.size[0] / 2, self.size[1] / 2, self.size[2] / 2)
+
+        self.scenario = 1  # 1 rising bubble, 2 RTI
+        self.config_dict = self.config()
+
+        self.csv_file = "benchmark.csv"
+
+    @wlb.member_callback
+    def config(self):
+        return {
+            'DomainSetup': {
+                'blocks': self.blocks,
+                'cellsPerBlock': self.cells,
+                'periodic': self.periodic,
+            },
+            'Parameters': {
+                'timesteps': self.timesteps,
+                'vtkWriteFrequency': self.vtkWriteFrequency,
+                'useGui': 0,
+                'remainingTimeLoggerFrequency': -1,
+                'timeStepStrategy': self.timeStepStrategy,
+                'overlappingWidth': self.overlappingWidth,
+                'gpuBlockSize': self.cuda_block_size,
+                'warmupSteps': self.warmupSteps,
+                'scenario': self.scenario,
+            },
+            'Boundaries_GPU': {
+                'Border': []
+            },
+            'Boundaries_CPU': {
+                'Border': []
+            },
+            'Bubble': {
+                'bubbleMidPoint': self.bubbleMidPoint,
+                'bubbleRadius': self.bubbleRadius,
+            },
+        }
+
+    @wlb.member_callback
+    def results_callback(self, **kwargs):
+        data = {}
+        data.update(self.config_dict['Parameters'])
+        data.update(self.config_dict['DomainSetup'])
+        data.update(kwargs)
+        data['executable'] = sys.argv[0]
+        data['compile_flags'] = wlb.build_info.compiler_flags
+        data['walberla_version'] = wlb.build_info.version
+        data['build_machine'] = wlb.build_info.build_machine
+        sequenceValuesToScalars(data)
+
+        df = pd.DataFrame.from_records([data])
+        if not os.path.isfile(self.csv_file):
+            df.to_csv(self.csv_file, index=False)
+        else:
+            df.to_csv(self.csv_file, index=False, mode='a', header=False)
+
+
+def overlap_benchmark():
+    scenarios = wlb.ScenarioManager()
+    overlappingWidths = [(1, 1, 1), (4, 1, 1), (8, 1, 1), (16, 1, 1), (32, 1, 1),
+                         (4, 4, 1), (8, 8, 1), (16, 16, 1), (32, 32, 1),
+                         (4, 4, 4), (8, 8, 8), (16, 16, 16), (32, 32, 32)]
+
+    cuda_blocks = [(32, 1, 1), (64, 1, 1), (128, 1, 1), (256, 1, 1), (512, 1, 1),
+                   (32, 2, 1), (64, 2, 1), (128, 2, 1), (256, 2, 1),
+                   (32, 4, 1), (64, 4, 1), (128, 4, 1),
+                   (32, 4, 2), (64, 4, 2), (128, 2, 2),
+                   (32, 8, 1), (64, 8, 1), (64, 4, 2),
+                   (32, 16, 1), (16, 16, 1), (16, 16, 2)]
+
+    # no overlap
+    scenarios.add(Scenario(timeStepStrategy='normal', overlappingWidth=(1, 1, 1), cuda_block_size=(16, 16, 1)))
+
+    # overlap
+    for overlap_strategy in ['overlap']:
+        for overlappingWidth in overlappingWidths:
+            for cuda_block in cuda_blocks:
+                scenario = Scenario(timeStepStrategy=overlap_strategy,
+                                    overlappingWidth=overlappingWidth,
+                                    cuda_block_size=cuda_block)
+                scenarios.add(scenario)
+
+
+def kernel_benchmark():
+    scenarios = wlb.ScenarioManager()
+
+    # overlap
+    # for overlap_strategy in ['phase_only', 'hydro_only', 'kernel_only']:
+    for overlap_strategy in ['overlap']:
+        scenario = Scenario(timeStepStrategy=overlap_strategy,
+                            overlappingWidth=(1, 1, 1),
+                            cuda_block_size=(128, 1, 1))
+        scenarios.add(scenario)
+
+
+# overlap_benchmark()
+kernel_benchmark()

apps/benchmarks/PhaseFieldAllenCahn/benchmark_cpu.py

+import os
+import waLBerla as wlb
+import pandas as pd
+
+from waLBerla.tools.sqlitedb import sequenceValuesToScalars
+
+import sys
+
+
+class Scenario:
+    def __init__(self, timeStepStrategy, overlappingWidth):
+        # output frequencies
+        self.vtkWriteFrequency = -1
+
+        # simulation parameters
+        self.timesteps = 500
+        edge_size = 50
+        self.cells = (edge_size, edge_size, edge_size)
+        self.blocks = (1, 1, 1)
+        self.periodic = (1, 1, 1)
+        self.size = (self.cells[0] * self.blocks[0],
+                     self.cells[1] * self.blocks[1],
+                     self.cells[2] * self.blocks[2])
+
+        self.timeStepStrategy = timeStepStrategy
+        self.overlappingWidth = overlappingWidth
+        self.cuda_block_size = (-1, -1, -1)
+        self.warmupSteps = 10
+
+        # bubble parameters
+        self.bubbleRadius = min(self.size) // 4
+        self.bubbleMidPoint = (self.size[0] / 2, self.size[1] / 2, self.size[2] / 2)
+
+        self.scenario = 1  # 1 rising bubble, 2 RTI
+        self.config_dict = self.config()
+
+        self.csv_file = "benchmark_cpu.csv"
+
+    @wlb.member_callback
+    def config(self):
+        return {
+            'DomainSetup': {
+                'blocks': self.blocks,
+                'cellsPerBlock': self.cells,
+                'periodic': self.periodic,
+            },
+            'Parameters': {
+                'timesteps': self.timesteps,
+                'vtkWriteFrequency': self.vtkWriteFrequency,
+                'useGui': 0,
+                'remainingTimeLoggerFrequency': 10.0,
+                'timeStepStrategy': self.timeStepStrategy,
+                'overlappingWidth': self.overlappingWidth,
+                'gpuBlockSize': self.cuda_block_size,
+                'warmupSteps': self.warmupSteps,
+                'scenario': self.scenario,
+            },
+            'Boundaries_GPU': {
+                'Border': []
+            },
+            'Boundaries_CPU': {
+                'Border': []
+            },
+            'Bubble': {
+                'bubbleMidPoint': self.bubbleMidPoint,
+                'bubbleRadius': self.bubbleRadius,
+            },
+        }
+
+    @wlb.member_callback
+    def results_callback(self, **kwargs):
+        data = {}
+        data.update(self.config_dict['Parameters'])
+        data.update(self.config_dict['DomainSetup'])
+        data.update(kwargs)
+        data['executable'] = sys.argv[0]
+        data['compile_flags'] = wlb.build_info.compiler_flags
+        data['walberla_version'] = wlb.build_info.version
+        data['build_machine'] = wlb.build_info.build_machine
+        sequenceValuesToScalars(data)
+
+        df = pd.DataFrame.from_records([data])
+        if not os.path.isfile(self.csv_file):
+            df.to_csv(self.csv_file, index=False)
+        else:
+            df.to_csv(self.csv_file, index=False, mode='a', header=False)
+
+
+def overlap_benchmark():
+    scenarios = wlb.ScenarioManager()
+    overlappingWidths = [(1, 1, 1), (4, 1, 1), (8, 1, 1), (16, 1, 1), (32, 1, 1),
+                         (4, 4, 1), (8, 8, 1), (16, 16, 1), (32, 32, 1),
+                         (4, 4, 4), (8, 8, 8), (16, 16, 16), (32, 32, 32),
+                         (64, 32, 32), (64, 64, 32), (64, 64, 64)]
+    # no overlap
+    scenarios.add(Scenario(timeStepStrategy='normal', overlappingWidth=(1, 1, 1)))
+
+    # overlap
+    for overlap_strategy in ['overlap']:
+        for overlappingWidth in overlappingWidths:
+            scenario = Scenario(timeStepStrategy=overlap_strategy,
+                                overlappingWidth=overlappingWidth)
+            scenarios.add(scenario)
+
+
+def kernel_benchmark():
+    scenarios = wlb.ScenarioManager()
+
+    # overlap
+    scenario = Scenario(timeStepStrategy='overlap',
+                        overlappingWidth=(8, 8, 8))
+    scenarios.add(scenario)
+
+
+# overlap_benchmark()
+kernel_benchmark()

apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.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 benchmark_multiphase.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformDirectScheme.h"
+
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+#include "core/math/Constants.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "field/AddToStorage.h"
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+#include "field/python/Exports.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "geometry/InitBoundaryHandling.h"
+
+#include "python_coupling/CreateConfig.h"
+#include "python_coupling/DictWrapper.h"
+#include "python_coupling/PythonCallback.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include <blockforest/communication/UniformBufferedScheme.h>
+
+#include "InitializerFunctions.h"
+
+//////////////////////////////
+// INCLUDE GENERATED FILES //
+////////////////////////////
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+#   include "cuda/AddGPUFieldToStorage.h"
+#   include "cuda/DeviceSelectMPI.h"
+#   include "cuda/HostFieldAllocator.h"
+#   include "cuda/NVTX.h"
+#   include "cuda/ParallelStreams.h"
+#   include "cuda/communication/GPUPackInfo.h"
+#   include "cuda/communication/MemcpyPackInfo.h"
+#   include "cuda/communication/UniformGPUScheme.h"
+
+#   include "GenDefines.h"
+#   include "PackInfo_phase_field.h"
+#   include "PackInfo_phase_field_distributions.h"
+#   include "PackInfo_velocity_based_distributions.h"
+#   include "hydro_LB_step.h"
+#   include "initialize_phase_field_distributions.h"
+#   include "initialize_velocity_based_distributions.h"
+#   include "phase_field_LB_step.h"
+
+#else
+#   include "GenDefines.h"
+#   include "PackInfo_phase_field.h"
+#   include "PackInfo_phase_field_distributions.h"
+#   include "PackInfo_velocity_based_distributions.h"
+#   include "hydro_LB_step.h"
+#   include "initialize_phase_field_distributions.h"
+#   include "initialize_velocity_based_distributions.h"
+#   include "phase_field_LB_step.h"
+#endif
+
+using namespace walberla;
+
+using PdfField_phase_T = GhostLayerField< real_t, Stencil_phase_T::Size >;
+using PdfField_hydro_T = GhostLayerField< real_t, Stencil_hydro_T::Size >;
+using VelocityField_T  = GhostLayerField< real_t, Stencil_hydro_T::Dimension >;
+using PhaseField_T     = GhostLayerField< real_t, 1 >;
+using flag_t           = walberla::uint8_t;
+using FlagField_T      = FlagField< flag_t >;
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+typedef cuda::GPUField< double > GPUField;
+#endif
+// using CommScheme_T = cuda::communication::UniformGPUScheme<stencil::D2Q9>;
+
+int main(int argc, char** argv)
+{
+   mpi::Environment env(argc, argv);
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+   cuda::selectDeviceBasedOnMpiRank();
+#endif
+
+   for (auto cfg = python_coupling::configBegin(argc, argv); cfg != python_coupling::configEnd(); ++cfg)
+   {
+      WALBERLA_MPI_WORLD_BARRIER()
+
+      auto config = *cfg;
+      logging::configureLogging(config);
+      shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGridFromConfig(config);
+
+      Vector3< uint_t > cellsPerBlock =
+         config->getBlock("DomainSetup").getParameter< Vector3< uint_t > >("cellsPerBlock");
+      // Reading parameters
+      auto parameters                    = config->getOneBlock("Parameters");
+      const std::string timeStepStrategy = parameters.getParameter< std::string >("timeStepStrategy", "normal");
+      const uint_t timesteps             = parameters.getParameter< uint_t >("timesteps", uint_c(50));
+      const real_t remainingTimeLoggerFrequency =
+         parameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0);
+      const uint_t scenario = parameters.getParameter< uint_t >("scenario", uint_c(1));
+      Vector3< int > overlappingWidth =
+         parameters.getParameter< Vector3< int > >("overlappingWidth", Vector3< int >(1, 1, 1));
+      const uint_t warmupSteps = parameters.getParameter< uint_t >("warmupSteps", uint_t(2));
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      // CPU fields
+      BlockDataID vel_field   = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
+      BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_t(0), field::fzyx);
+      // GPU fields
+      BlockDataID lb_phase_field_gpu = cuda::addGPUFieldToStorage< cuda::GPUField< real_t > >(
+         blocks, "lb phase field on GPU", Stencil_phase_T::Size, field::fzyx, 1);
+      BlockDataID lb_velocity_field_gpu = cuda::addGPUFieldToStorage< cuda::GPUField< real_t > >(
+         blocks, "lb velocity field on GPU", Stencil_hydro_T::Size, field::fzyx, 1);
+      BlockDataID vel_field_gpu =
+         cuda::addGPUFieldToStorage< VelocityField_T >(blocks, vel_field, "velocity field on GPU", true);
+      BlockDataID phase_field_gpu =
+         cuda::addGPUFieldToStorage< PhaseField_T >(blocks, phase_field, "phase field on GPU", true);
+#else
+      BlockDataID lb_phase_field =
+         field::addToStorage< PdfField_phase_T >(blocks, "lb phase field", real_t(0), field::fzyx);
+      BlockDataID lb_velocity_field =
+         field::addToStorage< PdfField_hydro_T >(blocks, "lb velocity field", real_t(0), field::fzyx);
+      BlockDataID vel_field   = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
+      BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_t(0), field::fzyx);
+#endif
+
+      if (timeStepStrategy != "phase_only" && timeStepStrategy != "hydro_only" && timeStepStrategy != "kernel_only")
+      {
+         WALBERLA_LOG_INFO_ON_ROOT("initialization of the phase field")
+         if (scenario == 1)
+         {
+            auto bubbleParameters = config->getOneBlock("Bubble");
+            const Vector3< real_t > bubbleMidPoint =
+               bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint");
+            const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", 20.0);
+            initPhaseField_bubble(blocks, phase_field, bubbleRadius, bubbleMidPoint);
+         }
+         else if (scenario == 2)
+         {
+            initPhaseField_RTI(blocks, phase_field);
+         }
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+         cuda::fieldCpy< GPUField, PhaseField_T >(blocks, phase_field_gpu, phase_field);
+#endif
+         WALBERLA_LOG_INFO_ON_ROOT("initialization of the phase field done")
+      }
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      Vector3< int32_t > gpuBlockSize =
+         parameters.getParameter< Vector3< int32_t > >("gpuBlockSize", Vector3< int32_t >(256, 1, 1));
+      pystencils::initialize_phase_field_distributions init_h(lb_phase_field_gpu, phase_field_gpu, vel_field_gpu);
+      pystencils::initialize_velocity_based_distributions init_g(lb_velocity_field_gpu, vel_field_gpu);
+
+      pystencils::phase_field_LB_step phase_field_LB_step(
+         lb_phase_field_gpu, phase_field_gpu, vel_field_gpu, gpuBlockSize[0], gpuBlockSize[1],
+         Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
+      pystencils::hydro_LB_step hydro_LB_step(lb_velocity_field_gpu, phase_field_gpu, vel_field_gpu, gpuBlockSize[0],
+                                              gpuBlockSize[1],
+                                              Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
+#else
+      pystencils::initialize_phase_field_distributions init_h(lb_phase_field, phase_field, vel_field);
+      pystencils::initialize_velocity_based_distributions init_g(lb_velocity_field, vel_field);
+      pystencils::phase_field_LB_step phase_field_LB_step(
+         lb_phase_field, phase_field, vel_field, Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
+      pystencils::hydro_LB_step hydro_LB_step(lb_velocity_field, phase_field, vel_field,
+                                              Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
+#endif
+
+// add communication
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      auto Comm_velocity_based_distributions =
+         make_shared< cuda::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, 0);
+      auto generatedPackInfo_velocity_based_distributions =
+         make_shared< pystencils::PackInfo_velocity_based_distributions >(lb_velocity_field_gpu);
+      Comm_velocity_based_distributions->addPackInfo(generatedPackInfo_velocity_based_distributions);
+      auto generatedPackInfo_phase_field = make_shared< pystencils::PackInfo_phase_field >(phase_field_gpu);
+      Comm_velocity_based_distributions->addPackInfo(generatedPackInfo_phase_field);
+
+      auto Comm_phase_field_distributions =
+         make_shared< cuda::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, 0);
+      auto generatedPackInfo_phase_field_distributions =
+         make_shared< pystencils::PackInfo_phase_field_distributions >(lb_phase_field_gpu);
+      Comm_phase_field_distributions->addPackInfo(generatedPackInfo_phase_field_distributions);
+#else
+
+
+      blockforest::communication::UniformBufferedScheme< Stencil_hydro_T > Comm_velocity_based_distributions(blocks);
+
+      auto generatedPackInfo_phase_field = make_shared< pystencils::PackInfo_phase_field >(phase_field);
+      auto generatedPackInfo_velocity_based_distributions =
+         make_shared< pystencils::PackInfo_velocity_based_distributions >(lb_velocity_field);
+
+      Comm_velocity_based_distributions.addPackInfo(generatedPackInfo_phase_field);
+      Comm_velocity_based_distributions.addPackInfo(generatedPackInfo_velocity_based_distributions);
+
+      blockforest::communication::UniformBufferedScheme< Stencil_hydro_T > Comm_phase_field_distributions(blocks);
+      auto generatedPackInfo_phase_field_distributions =
+         make_shared< pystencils::PackInfo_phase_field_distributions >(lb_phase_field);
+      Comm_phase_field_distributions.addPackInfo(generatedPackInfo_phase_field_distributions);
+#endif
+
+      BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+      // Boundaries
+      const FlagUID fluidFlagUID("Fluid");
+      auto boundariesConfig = config->getBlock("Boundaries_GPU");
+      if (boundariesConfig)
+      {
+         geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldID, boundariesConfig);
+         geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldID, fluidFlagUID);
+      }
+
+      // initialize the two lattice Boltzmann fields
+      if (timeStepStrategy != "phase_only" && timeStepStrategy != "hydro_only" && timeStepStrategy != "kernel_only")
+      {
+         WALBERLA_LOG_INFO_ON_ROOT("initialization of the distributions")
+         for (auto& block : *blocks)
+         {
+            init_h(&block);
+            init_g(&block);
+         }
+         WALBERLA_LOG_INFO_ON_ROOT("initialization of the distributions done")
+      }
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      int streamLowPriority  = 0;
+      int streamHighPriority = 0;
+      auto defaultStream     = cuda::StreamRAII::newPriorityStream(streamLowPriority);
+      auto innerOuterStreams = cuda::ParallelStreams(streamHighPriority);
+#endif
+
+      auto timeLoop = make_shared< SweepTimeloop >(blocks->getBlockStorage(), timesteps);
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      auto normalTimeStep = [&]() {
+         for (auto& block : *blocks)
+         {
+            Comm_phase_field_distributions->communicate(nullptr);
+            phase_field_LB_step(&block);
+
+            Comm_velocity_based_distributions->communicate(nullptr);
+            hydro_LB_step(&block);
+         }
+      };
+      auto simpleOverlapTimeStep = [&]() {
+         Comm_phase_field_distributions->startCommunication(defaultStream);
+         for (auto& block : *blocks)
+            phase_field_LB_step.inner(&block, defaultStream);
+         Comm_phase_field_distributions->wait(defaultStream);
+         for (auto& block : *blocks)
+            phase_field_LB_step.outer(&block, defaultStream);
+
+         Comm_velocity_based_distributions->startCommunication(defaultStream);
+         for (auto& block : *blocks)
+            hydro_LB_step.inner(&block, defaultStream);
+         Comm_velocity_based_distributions->wait(defaultStream);
+         for (auto& block : *blocks)
+            hydro_LB_step.outer(&block, defaultStream);
+      };
+      auto phase_only = [&]() {
+         for (auto& block : *blocks)
+            phase_field_LB_step(&block);
+      };
+      auto hydro_only = [&]() {
+         for (auto& block : *blocks)
+            hydro_LB_step(&block);
+      };
+      auto without_comm = [&]() {
+         for (auto& block : *blocks)
+            phase_field_LB_step(&block);
+         for (auto& block : *blocks)
+            hydro_LB_step(&block);
+      };
+#else
+      auto normalTimeStep = [&]() {
+         for (auto& block : *blocks)
+         {
+            Comm_phase_field_distributions();
+            phase_field_LB_step(&block);
+
+            Comm_velocity_based_distributions();
+            hydro_LB_step(&block);
+         }
+      };
+      auto simpleOverlapTimeStep = [&]() {
+        Comm_phase_field_distributions.startCommunication();
+        for (auto& block : *blocks)
+           phase_field_LB_step.inner(&block);
+        Comm_phase_field_distributions.wait();
+        for (auto& block : *blocks)
+           phase_field_LB_step.outer(&block);
+
+        Comm_velocity_based_distributions.startCommunication();
+        for (auto& block : *blocks)
+           hydro_LB_step.inner(&block);
+        Comm_velocity_based_distributions.wait();
+        for (auto& block : *blocks)
+           hydro_LB_step.outer(&block);
+      };
+      auto phase_only = [&]() {
+         for (auto& block : *blocks)
+            phase_field_LB_step(&block);
+      };
+      auto hydro_only = [&]() {
+         for (auto& block : *blocks)
+            hydro_LB_step(&block);
+      };
+      auto without_comm = [&]() {
+         for (auto& block : *blocks)
+            phase_field_LB_step(&block);
+         for (auto& block : *blocks)
+            hydro_LB_step(&block);
+      };
+#endif
+      std::function< void() > timeStep;
+      if (timeStepStrategy == "overlap")
+      {
+         timeStep = std::function< void() >(simpleOverlapTimeStep);
+         WALBERLA_LOG_INFO_ON_ROOT("overlapping timestep")
+      }
+      else if (timeStepStrategy == "phase_only")
+      {
+         timeStep = std::function< void() >(phase_only);
+         WALBERLA_LOG_INFO_ON_ROOT("started only phasefield step without communication for benchmarking")
+      }
+      else if (timeStepStrategy == "hydro_only")
+      {
+         timeStep = std::function< void() >(hydro_only);
+         WALBERLA_LOG_INFO_ON_ROOT("started only hydro step without communication for benchmarking")
+      }
+      else if (timeStepStrategy == "kernel_only")
+      {
+         timeStep = std::function< void() >(without_comm);
+         WALBERLA_LOG_INFO_ON_ROOT("started complete phasefield model without communication for benchmarking")
+      }
+      else
+      {
+         timeStep = std::function< void() >(normalTimeStep);
+         WALBERLA_LOG_INFO_ON_ROOT("normal timestep with no overlapping")
+      }
+
+      timeLoop->add() << BeforeFunction(timeStep) << Sweep([](IBlock*) {}, "time step");
+
+      // remaining time logger
+      if (remainingTimeLoggerFrequency > 0)
+         timeLoop->addFuncAfterTimeStep(
+            timing::RemainingTimeLogger(timeLoop->getNrOfTimeSteps(), remainingTimeLoggerFrequency),
+            "remaining time logger");
+
+      uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0);
+      if (vtkWriteFrequency > 1)
+      {
+         auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false, "vtk_out",
+                                                         "simulation_step", false, true, true, false, 0);
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+         vtkOutput->addBeforeFunction([&]() {
+            cuda::fieldCpy< PhaseField_T, GPUField >(blocks, phase_field, phase_field_gpu);
+            // cuda::fieldCpy<VelocityField_T, GPUField>( blocks, vel_field, vel_field_gpu );
+         });
+#endif
+         auto phaseWriter = make_shared< field::VTKWriter< PhaseField_T > >(phase_field, "phase");
+         vtkOutput->addCellDataWriter(phaseWriter);
+
+         timeLoop->addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
+      }
+
+      for (uint_t i = 0; i < warmupSteps; ++i)
+         timeLoop->singleStep();
+
+      timeLoop->setCurrentTimeStepToZero();
+      WALBERLA_LOG_INFO_ON_ROOT("Starting simulation with " << timesteps << " time steps")
+      WcTimer simTimer;
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      cudaDeviceSynchronize();
+#endif
+      simTimer.start();
+      timeLoop->run();
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      cudaDeviceSynchronize();
+#endif
+      simTimer.end();
+      WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
+      auto time            = simTimer.last();
+      auto nrOfCells       = real_c(cellsPerBlock[0] * cellsPerBlock[1] * cellsPerBlock[2]);
+      auto mlupsPerProcess = nrOfCells * real_c(timesteps) / time * 1e-6;
+      WALBERLA_LOG_RESULT_ON_ROOT("MLUPS per process: " << mlupsPerProcess)
+      WALBERLA_LOG_RESULT_ON_ROOT("Time per time step: " << time / real_c(timesteps) << " s")
+      WALBERLA_ROOT_SECTION()
+      {
+         python_coupling::PythonCallback pythonCallbackResults("results_callback");
+         if (pythonCallbackResults.isCallable())
+         {
+            pythonCallbackResults.data().exposeValue("mlupsPerProcess", mlupsPerProcess);
+            // Call Python function to report results
+            pythonCallbackResults();
+         }
+      }
+   }
+
+   return EXIT_SUCCESS;
+}

apps/benchmarks/PhaseFieldAllenCahn/benchmark_piz_daint.py

+import os
+import waLBerla as wlb
+import csv
+
+from waLBerla.tools.sqlitedb import sequenceValuesToScalars
+
+import sys
+
+
+class Scenario:
+    def __init__(self, timeStepStrategy, overlappingWidth, cuda_block_size):
+        # output frequencies
+        self.vtkWriteFrequency = 0
+
+        # simulation parameters
+        self.timesteps = 201
+        edge_size = 100
+        self.cells = (edge_size, edge_size, edge_size)
+        self.blocks = (1, 1, 1)
+        self.periodic = (1, 1, 1)
+        self.size = (self.cells[0] * self.blocks[0],
+                     self.cells[1] * self.blocks[1],
+                     self.cells[2] * self.blocks[2])
+
+        self.timeStepStrategy = timeStepStrategy
+        self.overlappingWidth = overlappingWidth
+        self.cuda_block_size = cuda_block_size
+        self.warmupSteps = 20
+
+        # bubble parameters
+        self.bubbleRadius = min(self.size) // 4
+        self.bubbleMidPoint = (self.size[0] / 2, self.size[1] / 2, self.size[2] / 2)
+
+        self.scenario = 1  # 1 rising bubble, 2 RTI
+        self.config_dict = self.config()
+
+        self.csv_file = "benchmark_piz_daint.csv"
+
+    @wlb.member_callback
+    def config(self):
+        return {
+            'DomainSetup': {
+                'blocks': self.blocks,
+                'cellsPerBlock': self.cells,
+                'periodic': self.periodic,
+            },
+            'Parameters': {
+                'timesteps': self.timesteps,
+                'vtkWriteFrequency': self.vtkWriteFrequency,
+                'useGui': 0,
+                'remainingTimeLoggerFrequency': -1,
+                'timeStepStrategy': self.timeStepStrategy,
+                'overlappingWidth': self.overlappingWidth,
+                'gpuBlockSize': self.cuda_block_size,
+                'warmupSteps': self.warmupSteps,
+                'scenario': self.scenario,
+            },
+            'Boundaries_GPU': {
+                'Border': []
+            },
+            'Boundaries_CPU': {
+                'Border': []
+            },
+            'Bubble': {
+                'bubbleMidPoint': self.bubbleMidPoint,
+                'bubbleRadius': self.bubbleRadius,
+            },
+        }
+
+    @wlb.member_callback
+    def results_callback(self, **kwargs):
+        data = {}
+        data.update(self.config_dict['Parameters'])
+        data.update(self.config_dict['DomainSetup'])
+        data.update(kwargs)
+        data['executable'] = sys.argv[0]
+        data['compile_flags'] = wlb.build_info.compiler_flags
+        data['walberla_version'] = wlb.build_info.version
+        data['build_machine'] = wlb.build_info.build_machine
+        sequenceValuesToScalars(data)
+
+        if not os.path.isfile(self.csv_file):
+            try:
+                with open(self.csv_file, 'w') as csvfile:
+                    writer = csv.DictWriter(csvfile, fieldnames=data)
+                    writer.writeheader()
+                    writer.writerow(data)
+            except IOError:
+                print("could not create csv file")
+        else:
+            try:
+                with open(self.csv_file, 'a') as csvfile:
+                    writer = csv.DictWriter(csvfile, fieldnames=data)
+                    writer.writerow(data)
+            except IOError:
+                print("could not write to csv file")
+
+
+def overlap_benchmark():
+    scenarios = wlb.ScenarioManager()
+    overlappingWidths = [(1, 1, 1), (4, 1, 1), (8, 1, 1), (16, 1, 1), (32, 1, 1),
+                         (4, 4, 1), (8, 8, 1), (16, 16, 1), (32, 32, 1),
+                         (4, 4, 4), (8, 8, 8), (16, 16, 16), (32, 32, 32)]
+
+    cuda_blocks = [(32, 1, 1), (64, 1, 1), (128, 1, 1), (256, 1, 1),
+                   (32, 2, 1), (64, 2, 1), (128, 2, 1),
+                   (32, 4, 1), (64, 4, 1),
+                   (32, 4, 2), (32, 8, 1), (16, 16, 1)]
+
+    # no overlap
+    scenarios.add(Scenario(timeStepStrategy='normal', overlappingWidth=(1, 1, 1), cuda_block_size=(16, 16, 1)))
+
+    # overlap
+    for overlap_strategy in ['overlap']:
+        for overlappingWidth in overlappingWidths:
+            for cuda_block in cuda_blocks:
+                scenario = Scenario(timeStepStrategy=overlap_strategy,
+                                    overlappingWidth=overlappingWidth,
+                                    cuda_block_size=cuda_block)
+                scenarios.add(scenario)
+
+
+def kernel_benchmark():
+    scenarios = wlb.ScenarioManager()
+
+    # overlap
+    # for overlap_strategy in ['phase_only', 'hydro_only', 'kernel_only']:
+    for overlap_strategy in ['overlap']:
+        scenario = Scenario(timeStepStrategy=overlap_strategy,
+                            overlappingWidth=(16, 16, 16),
+                            cuda_block_size=(128, 4, 1))
+        scenarios.add(scenario)
+
+
+def weak_scaling():
+    scenarios = wlb.ScenarioManager()
+
+    # overlap
+    # for overlap_strategy in ['phase_only', 'hydro_only', 'kernel_only']:
+    for overlap_strategy in ['overlap']:
+        scenario = Scenario(timeStepStrategy=overlap_strategy,
+                            overlappingWidth=(32, 32, 32),
+                            cuda_block_size=(128, 2, 1))
+        scenarios.add(scenario)
+
+
+# overlap_benchmark()
+# kernel_benchmark()
+weak_scaling()

apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py

+from pystencils import fields, TypedSymbol
+from pystencils.simp import sympy_cse
+from pystencils import AssignmentCollection
+
+from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
+from lbmpy.stencils import get_stencil
+
+from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_from_kernel
+
+from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_phase_field_lb, \
+    initializer_kernel_hydro_lb, interface_tracking_force, \
+    hydrodynamic_force, get_collision_assignments_hydro
+
+from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
+
+import numpy as np
+
+stencil_phase = get_stencil("D3Q15", "walberla")
+stencil_hydro = get_stencil("D3Q27", "walberla")
+q_phase = len(stencil_phase)
+q_hydro = len(stencil_hydro)
+
+maxwellian_moments = False
+assert (len(stencil_phase[0]) == len(stencil_hydro[0]))
+dimensions = len(stencil_hydro[0])
+
+########################
+# PARAMETER DEFINITION #
+########################
+
+density_liquid = 1.0
+density_gas = 0.001
+surface_tension = 0.0001
+M = 0.02
+
+# phase-field parameter
+drho3 = (density_liquid - density_gas) / 3
+# interface thickness
+W = 5
+# coefficient related to surface tension
+beta = 12.0 * (surface_tension / W)
+# coefficient related to surface tension
+kappa = 1.5 * surface_tension * W
+# relaxation rate allen cahn (h)
+w_c = 1.0 / (0.5 + (3.0 * M))
+
+########################
+# FIELDS #
+########################
+
+# velocity field
+u = fields(f"vel_field({dimensions}): [{dimensions}D]", layout='fzyx')
+# phase-field
+C = fields(f"phase_field: [{dimensions}D]", layout='fzyx')
+
+# phase-field distribution functions
+h = fields(f"lb_phase_field({q_phase}): [{dimensions}D]", layout='fzyx')
+h_tmp = fields(f"lb_phase_field_tmp({q_phase}): [{dimensions}D]", layout='fzyx')
+# hydrodynamic distribution functions
+g = fields(f"lb_velocity_field({q_hydro}): [{dimensions}D]", layout='fzyx')
+g_tmp = fields(f"lb_velocity_field_tmp({q_hydro}): [{dimensions}D]", layout='fzyx')
+
+########################################
+# RELAXATION RATES AND EXTERNAL FORCES #
+########################################
+
+# calculate the relaxation rate for the hydro lb as well as the body force
+density = density_gas + C.center * (density_liquid - density_gas)
+
+# force acting on all phases of the model
+body_force = np.array([0, 1e-6, 0])
+
+relaxation_time = 0.03 + 0.5
+relaxation_rate = 1.0 / relaxation_time
+
+###############
+# LBM METHODS #
+###############
+
+method_phase = create_lb_method(stencil=stencil_phase, method='srt', relaxation_rate=w_c, compressible=True)
+
+method_hydro = create_lb_method(stencil=stencil_hydro, method="mrt", weighted=True,
+                                relaxation_rates=[relaxation_rate, 1, 1, 1, 1, 1],
+                                maxwellian_moments=maxwellian_moments)
+
+# create the kernels for the initialization of the g and h field
+h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
+g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
+
+
+force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
+force_model_h = MultiphaseForceModel(force=force_h)
+
+force_g = hydrodynamic_force(g, C, method_hydro, relaxation_time, density_liquid, density_gas, kappa, beta, body_force)
+
+h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
+sum_h = np.sum(h_tmp_symbol_list[:])
+
+####################
+# LBM UPDATE RULES #
+####################
+
+method_phase.set_force_model(force_model_h)
+
+phase_field_LB_step = create_lb_update_rule(lb_method=method_phase,
+                                            velocity_input=u,
+                                            compressible=True,
+                                            optimization={"symbolic_field": h,
+                                                          "symbolic_temporary_field": h_tmp},
+                                            kernel_type='stream_pull_collide')
+
+
+phase_field_LB_step.set_main_assignments_from_dict({**phase_field_LB_step.main_assignments_dict, **{C.center: sum_h}})
+
+phase_field_LB_step = AssignmentCollection(main_assignments=phase_field_LB_step.main_assignments,
+                                           subexpressions=phase_field_LB_step.subexpressions)
+phase_field_LB_step = sympy_cse(phase_field_LB_step)
+
+# ---------------------------------------------------------------------------------------------------------
+
+hydro_LB_step = get_collision_assignments_hydro(lb_method=method_hydro,
+                                                density=density,
+                                                velocity_input=u,
+                                                force=force_g,
+                                                optimization={"symbolic_field": g,
+                                                              "symbolic_temporary_field": g_tmp},
+                                                kernel_type='collide_stream_push')
+
+# streaming of the hydrodynamic distribution
+stream_hydro = create_lb_update_rule(stencil=stencil_hydro,
+                                     optimization={"symbolic_field": g,
+                                                   "symbolic_temporary_field": g_tmp},
+                                     kernel_type='stream_pull_only')
+
+###################
+# GENERATE SWEEPS #
+###################
+
+cpu_vec = {'instruction_set': 'sse', 'assume_inner_stride_one': True, 'nontemporal': True}
+
+vp = [('int32_t', 'cudaBlockSize0'),
+      ('int32_t', 'cudaBlockSize1')]
+
+sweep_block_size = (TypedSymbol("cudaBlockSize0", np.int32),
+                    TypedSymbol("cudaBlockSize1", np.int32),
+                    1)
+sweep_params = {'block_size': sweep_block_size}
+
+info_header = f"""
+#include "stencil/D3Q{q_phase}.h"\nusing Stencil_phase_T = walberla::stencil::D3Q{q_phase}; 
+#include "stencil/D3Q{q_hydro}.h"\nusing Stencil_hydro_T = walberla::stencil::D3Q{q_hydro};
+"""
+
+with CodeGeneration() as ctx:
+    if ctx.cuda is False:
+        if not ctx.optimize_for_localhost:
+            cpu_vec = {'instruction_set': None}
+
+        generate_sweep(ctx, 'initialize_phase_field_distributions', h_updates)
+        generate_sweep(ctx, 'initialize_velocity_based_distributions', g_updates)
+
+        generate_sweep(ctx, 'phase_field_LB_step', phase_field_LB_step,
+                       field_swaps=[(h, h_tmp)],
+                       inner_outer_split=True,
+                       cpu_vectorize_info=cpu_vec)
+
+        generate_sweep(ctx, 'hydro_LB_step', hydro_LB_step,
+                       field_swaps=[(g, g_tmp)],
+                       inner_outer_split=True,
+                       cpu_vectorize_info=cpu_vec)
+
+        # communication
+        generate_pack_info_from_kernel(ctx, 'PackInfo_phase_field_distributions',
+                                       phase_field_LB_step.main_assignments, target='cpu')
+        generate_pack_info_from_kernel(ctx, 'PackInfo_phase_field',
+                                       hydro_LB_step.all_assignments, target='cpu', kind='pull')
+        generate_pack_info_from_kernel(ctx, 'PackInfo_velocity_based_distributions',
+                                       hydro_LB_step.all_assignments, target='cpu', kind='push')
+
+        ctx.write_file("GenDefines.h", info_header)
+
+    if ctx.cuda:
+        generate_sweep(ctx, 'initialize_phase_field_distributions',
+                       h_updates, target='gpu')
+        generate_sweep(ctx, 'initialize_velocity_based_distributions',
+                       g_updates, target='gpu')
+
+        generate_sweep(ctx, 'phase_field_LB_step', phase_field_LB_step,
+                       field_swaps=[(h, h_tmp)],
+                       inner_outer_split=True,
+                       target='gpu',
+                       gpu_indexing_params=sweep_params,
+                       varying_parameters=vp)
+
+        generate_sweep(ctx, 'hydro_LB_step', hydro_LB_step,
+                       field_swaps=[(g, g_tmp)],
+                       inner_outer_split=True,
+                       target='gpu',
+                       gpu_indexing_params=sweep_params,
+                       varying_parameters=vp)
+        # communication
+        generate_pack_info_from_kernel(ctx, 'PackInfo_phase_field_distributions',
+                                       phase_field_LB_step.main_assignments, target='gpu')
+        generate_pack_info_from_kernel(ctx, 'PackInfo_phase_field',
+                                       hydro_LB_step.all_assignments, target='gpu', kind='pull')
+        generate_pack_info_from_kernel(ctx, 'PackInfo_velocity_based_distributions',
+                                       hydro_LB_step.all_assignments, target='gpu', kind='push')
+
+        ctx.write_file("GenDefines.h", info_header)
+
+print("finished code generation successfully")

apps/benchmarks/PhaseFieldAllenCahn/profiling.py

+import os
+import waLBerla as wlb
+
+
+class Scenario:
+    def __init__(self):
+        # output frequencies
+        self.vtkWriteFrequency = 100
+
+        # simulation parameters
+        self.timesteps = 3
+        edge_size = 200
+        self.cells = (edge_size, edge_size, edge_size)
+        self.blocks = (1, 1, 1)
+        self.periodic = (1, 1, 1)
+        self.size = (self.cells[0] * self.blocks[0],
+                     self.cells[1] * self.blocks[1],
+                     self.cells[2] * self.blocks[2])
+
+        self.timeStepStrategy = 'phase_only'
+        self.overlappingWidth = (1, 1, 1)
+        self.cuda_block_size = (128, 1, 1)
+
+        # bubble parameters
+        self.bubbleRadius = min(self.size) // 4
+        self.bubbleMidPoint = (self.size[0] / 2, self.size[1] / 2, self.size[2] / 2)
+
+        self.scenario = 1  # 1 rising bubble, 2 RTI
+        self.config_dict = self.config()
+
+        self.csv_file = "benchmark.csv"
+
+    @wlb.member_callback
+    def config(self):
+        return {
+            'DomainSetup': {
+                'blocks': self.blocks,
+                'cellsPerBlock': self.cells,
+                'periodic': self.periodic,
+            },
+            'Parameters': {
+                'timesteps': self.timesteps,
+                'vtkWriteFrequency': self.vtkWriteFrequency,
+                'useGui': 0,
+                'remainingTimeLoggerFrequency': -1,
+                'timeStepStrategy': self.timeStepStrategy,
+                'overlappingWidth': self.overlappingWidth,
+                'gpuBlockSize': self.cuda_block_size,
+                'warmupSteps': 0,
+                'scenario': self.scenario,
+            },
+            'Boundaries_GPU': {
+                'Border': []
+            },
+            'Boundaries_CPU': {
+                'Border': []
+            },
+            'Bubble': {
+                'bubbleMidPoint': self.bubbleMidPoint,
+                'bubbleRadius': self.bubbleRadius,
+            },
+        }
+
+
+scenarios = wlb.ScenarioManager()
+scenarios.add(Scenario())

apps/showcases/CMakeLists.txt

@@ -2,3 +2,6 @@ add_subdirectory( BidisperseFluidizedBed )
 add_subdirectory( CombinedResolvedUnresolved )
 add_subdirectory( HeatConduction )
 add_subdirectory( Mixer )
+if ( WALBERLA_BUILD_WITH_CODEGEN)
+add_subdirectory( PhaseFieldAllenCahn )
+endif()

apps/showcases/PhaseFieldAllenCahn/CMakeLists.txt

+if ( NOT WALBERLA_BUILD_WITH_CUDA)
+add_subdirectory( CPU )
+endif()
+
+if ( WALBERLA_BUILD_WITH_CUDA)
+add_subdirectory( GPU )
+endif()

apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt

+waLBerla_link_files_to_builddir(*.prm)
+waLBerla_link_files_to_builddir(*.py)
+waLBerla_link_files_to_builddir(*.obj)
+
+waLBerla_generate_target_from_python(NAME PhaseFieldCodeGenCPU
+        FILE multiphase_codegen.py
+        OUT_FILES initialize_phase_field_distributions.cpp initialize_phase_field_distributions.h
+        initialize_velocity_based_distributions.cpp initialize_velocity_based_distributions.h
+        phase_field_LB_step.cpp phase_field_LB_step.h
+        phase_field_LB_NoSlip.cpp phase_field_LB_NoSlip.h
+        hydro_LB_step.cpp hydro_LB_step.h
+        hydro_LB_NoSlip.cpp hydro_LB_NoSlip.h
+        stream_hydro.cpp stream_hydro.h
+        GenDefines.h)
+
+waLBerla_add_executable(NAME multiphase
+        FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp contact.cpp CalculateNormals.cpp multiphase_codegen.py
+        DEPENDS blockforest core field postprocessing lbm geometry timeloop gui PhaseFieldCodeGenCPU)
+
+
+

apps/showcases/PhaseFieldAllenCahn/CPU/CalculateNormals.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 CalculateNormals.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "CalculateNormals.h"
+
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+
+#include "field/FlagField.h"
+
+namespace walberla
+{
+using FlagField_T    = FlagField< uint8_t >;
+using NormalsField_T = GhostLayerField< int8_t, 3 >;
+
+void calculate_normals(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID normalsFieldID,
+                       ConstBlockDataID flagFieldID, FlagUID domainFlagUID, FlagUID boundaryFlagUID)
+{
+   for (auto& block : *blocks)
+   {
+      CellInterval globalCellBB = blocks->getBlockCellBB(block);
+      CellInterval blockLocalCellBB;
+      blocks->transformGlobalToBlockLocalCellInterval(blockLocalCellBB, block, globalCellBB);
+
+      auto* normalsField = block.getData< NormalsField_T >(normalsFieldID);
+      auto* flagField    = block.getData< FlagField_T >(flagFieldID);
+      auto boundaryFlag  = flagField->getFlag(boundaryFlagUID);
+      auto domainFlag    = flagField->getFlag(domainFlagUID);
+
+      // clang-format off
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(normalsField,
+
+         if( x < blockLocalCellBB.xMax() ){
+            if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x + 1, y, z) == domainFlag)
+               normalsField->get(x, y, z, 0) = 1;
+         }
+
+         if( x > blockLocalCellBB.xMin() ){
+            if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x - 1, y, z) == domainFlag)
+               normalsField->get(x, y, z, 0) = - 1;
+         }
+
+         if( y < blockLocalCellBB.yMax() ){
+            if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y + 1, z) == domainFlag)
+               normalsField->get(x, y, z, 1) = 1;
+         }
+
+         if( y > blockLocalCellBB.yMin() ){
+            if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y - 1, z) == domainFlag)
+               normalsField->get(x, y, z, 1) = - 1;
+         }
+
+         if( z < blockLocalCellBB.zMax() ){
+            if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y, z + 1) == domainFlag)
+               normalsField->get(x, y, z, 2) = 1;
+         }
+
+         if( z > blockLocalCellBB.zMin() ){
+            if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y, z - 1) == domainFlag)
+               normalsField->get(x, y, z, 2) = - 1;
+         }
+
+      )
+      // clang-format on
+   }
+}
+} // namespace walberla

apps/showcases/PhaseFieldAllenCahn/CPU/CalculateNormals.h

+//======================================================================================================================
+//
+//  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 CalculateNormals.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/StructuredBlockForest.h"
+
+#include "domain_decomposition/BlockDataID.h"
+#include "domain_decomposition/IBlock.h"
+
+#include "field/FlagField.h"
+
+namespace walberla
+{
+void calculate_normals(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID normalsFieldID,
+                       ConstBlockDataID flagFieldID, FlagUID fluidFlagUID, FlagUID boundaryFlagUID);
+}
\ No newline at end of file

apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.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 InitializerFunctions.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "core/logging/Initialization.h"
+#include "core/math/Constants.h"
+#include "core/math/Random.h"
+
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "python_coupling/DictWrapper.h"
+
+namespace walberla
+{
+using PhaseField_T = GhostLayerField< real_t, 1 >;
+using FlagField_T  = FlagField< uint8_t >;
+
+void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
+                           const real_t R                         = 10,
+                           const Vector3< real_t > bubbleMidPoint = Vector3< real_t >(0.0, 0.0, 0.0),
+                           const bool bubble = true, const real_t W = 5)
+{
+   for (auto& block : *blocks)
+   {
+      auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
+      // clang-format off
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell;
+         blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+         real_t Ri = sqrt((globalCell[0] - bubbleMidPoint[0]) * (globalCell[0] - bubbleMidPoint[0]) +
+                          (globalCell[1] - bubbleMidPoint[1]) * (globalCell[1] - bubbleMidPoint[1]) +
+                          (globalCell[2] - bubbleMidPoint[2]) * (globalCell[2] - bubbleMidPoint[2]));
+         if (bubble) phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (Ri - R) / W);
+         else phaseField->get(x, y, z)        = 0.5 - 0.5 * tanh(2.0 * (Ri - R) / W);
+      )
+      // clang-format on
+   }
+}
+
+void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
+                        const real_t W = 5)
+{
+   auto X              = blocks->getDomainCellBB().xMax();
+   auto halfY          = (blocks->getDomainCellBB().yMax()) / 2.0;
+   double perturbation = 0.05;
+
+   for (auto& block : *blocks)
+   {
+      auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
+      // clang-format off
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell;
+         blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+         real_t tmp =
+         perturbation * X * (cos((2.0 * math::pi * globalCell[0]) / X) + cos((2.0 * math::pi * globalCell[2]) / X));
+         phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) - tmp) / (W / 2.0));
+      )
+      // clang-format on
+   }
+}
+} // namespace walberla

apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h

+//======================================================================================================================
+//
+//  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 InitializerFunctions.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+#include "core/math/Constants.h"
+
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "python_coupling/DictWrapper.h"
+#pragma once
+
+namespace walberla
+{
+void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
+                           Vector3< real_t > bubbleMidPoint, bool bubble = true, real_t W = 5);
+
+void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W = 5);
+
+} // namespace walberla

apps/showcases/PhaseFieldAllenCahn/CPU/PythonExports.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 PythonExports.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "waLBerlaDefinitions.h"
+#ifdef WALBERLA_BUILD_WITH_PYTHON
+
+#include "python_coupling/Manager.h"
+
+#include "core/math/Constants.h"
+
+#include "blockforest/python/Exports.h"
+
+#include "field/communication/PackInfo.h"
+#include "field/python/Exports.h"
+
+#include "geometry/python/Exports.h"
+
+#include "postprocessing/python/Exports.h"
+
+#include "timeloop/python/Exports.h"
+
+
+namespace walberla {
+    using flag_t = uint8_t;
+    // clang-format off
+    void exportDataStructuresToPython() {
+
+        namespace bmpl = boost::mpl;
+
+        auto pythonManager = python_coupling::Manager::instance();
+
+        typedef bmpl::vector<
+                Field<walberla::real_t, 1>,
+                Field<walberla::real_t, 3>,
+                Field<walberla::real_t, 9>,
+                Field<walberla::real_t, 19>,
+                Field<walberla::real_t, 27>>
+                FieldTypes;
+
+        typedef bmpl::vector<stencil::D2Q9,
+                stencil::D3Q19,
+                stencil::D3Q27>
+                Stencils;
+
+        typedef bmpl::vector<
+                GhostLayerField<real_t,1>,
+                GhostLayerField<real_t,3>>
+                RealFieldTypes;
+
+        typedef bmpl::vector<
+                FlagField<flag_t>>
+                FlagFieldTypes;
+        // Field
+        pythonManager->addExporterFunction(field::exportModuleToPython<FieldTypes>);
+        pythonManager->addExporterFunction(field::exportGatherFunctions<FieldTypes>);
+        pythonManager->addBlockDataConversion<FieldTypes>();
+
+        // Blockforest
+        pythonManager->addExporterFunction(blockforest::exportModuleToPython<Stencils>);
+
+        // Timeloop
+        pythonManager->addExporterFunction(timeloop::exportModuleToPython);
+
+        // Postprocessing
+        pythonManager->addExporterFunction( postprocessing::exportModuleToPython<RealFieldTypes, FlagFieldTypes> );
+
+        // Geometry
+        pythonManager->addExporterFunction( geometry::exportModuleToPython );
+    }
+   // clang-format on
+}
+
+#else
+
+namespace walberla {
+   void exportDataStructuresToPython() {}
+} // namespace walberla
+
+#endif
\ No newline at end of file

apps/showcases/PhaseFieldAllenCahn/CPU/PythonExports.h

+//======================================================================================================================
+//
+//  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 PythonExports.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+namespace walberla
+{
+void exportDataStructuresToPython();
+} // namespace walberla